Magnetic resonance spectroscopy pulse sequence, acquisition, and processing system and method

ABSTRACT

Systems and methods are provided for processing a set of multiple serially acquired magnetic resonance spectroscopy (MRS) free induction decay (FID) frames from a multi-frame MRS acquisition series from a region of interest (ROI) in a subject, and for providing a post-processed MRS spectrum. Processing parameters are dynamically varied while measuring results to determine the optimal post-processed results. Spectral regions opposite water from chemical regions of interest are evaluated and used in at least one processing operation. Frequency shift error is estimated via spectral correlation between free induction decay (FID) frames and a reference spectrum. Multiple groups of FID frames within the acquired set are identified to different phases corresponding with a phase step cycle of the acquisition. Baseline correction is also performed via rank order filter (ROF) estimate and a polynomial fit. Sections of the ROF may be excluded from the polynomial fit, such as for example sections determined to be associated with relevant spectral peaks.

CROSS-REFERENCE TO RELATED APPLICATIONS

This application is a continuation of U.S. patent application Ser. No.16/100,788, filed Aug. 10, 2018, and titled MAGNETIC RESONANCESPECTROSCOPY PULSE SEQUENCE, ACQUISITION, AND PROCESSING SYSTEM ANDMETHOD, which is a continuation of U.S. patent application Ser. No.15/160,423, filed May 20, 2016, and titled MAGNETIC RESONANCESPECTROSCOPY PULSE SEQUENCE, ACQUISITION, AND PROCESSING SYSTEM ANDMETHOD, which is a continuation of U.S. patent application Ser. No.14/625,918, filed Feb. 19, 2015, and titled MAGNETIC RESONANCESPECTROSCOPY PULSE SEQUENCE, ACQUISITION, AND PROCESSING SYSTEM ANDMETHOD, which is a continuation of U.S. patent application Ser. No.13/843,117, filed Mar. 15, 2013, and titled MAGNETIC RESONANCESPECTROSCOPY PULSE SEQUENCE, ACQUISITION, AND PROCESSING SYSTEM ANDMETHOD, which claims the benefit under 35 U.S.C. § 119(e) of U.S.Provisional Patent Application No. 61/624,284, filed Apr. 14, 2012, andtitled MAGNETIC RESONANCE SPECTROSCOPY PULSE SEQUENCE, ACQUISITION, ANDPROCESSING SYSTEM AND METHOD. The above-identified applications arehereby incorporated by reference in their entirety and made a part ofthis specification for all that they disclose.

INCORPORATION BY REFERENCE

The following disclosures are hereby incorporated by reference in theirentirety and made a part of this specification for all that theydisclose: U.S. Patent Publication No. 2008/0039710, filed Jul. 27, 2007,and titled “SYSTEM AND METHODS USING NUCLEAR MAGNETIC RESONANCE (NMR)SPECTROSCOPY TO EVALUATE PAIN AND DEGENERATIVE PROPERTIES OF TISSUE”;U.S. Patent Publication No. 2009/0030308, filed Mar. 21, 2008, andtitled “SYSTEM, COMPOSITION, AND METHODS FOR LOCAL IMAGING AND TREATMENTOF PAIN”; International Patent Publication No. WO 2009/148550, filed May29, 2009, and titled “BIOMARKERS FOR PAINFUL INTERVERTEBRAL DISCS ANDMETHODS OF USE THEREOF”; U.S. Patent Publication No. 2011/0087087, filedOct. 14, 2009, and titled “MR SPECTROSCOPYSYSTEM AND METHOD FORDIAGNOSING PAINFUL AND NON-PAINFUL INTERVERTEBRAL DISCS”; InternationalPatent Publication No. WO 2011/047197, filed Oct. 14, 2010, and titled“MR SPECTROSCOPY SYSTEM AND METHOD FOR DIAGNOSING PAINFUL ANDNON-PAINFUL INTERVERTEBRAL DISCS”; and U.S. patent application Ser. No.13/830,632, filed Mar. 14, 2013, and titled SYSTEMS AND METHODS FORAUTOMATED VOXELATION OF REGIONS OF INTEREST FOR MAGNETIC RESONANCESPECTROSCOPY.

BACKGROUND OF THE INVENTION Field of the Invention

This disclosure relates to systems, processors, devices, and methods formeasuring chemical constituents in tissue for diagnosing medicalconditions. More specifically, it relates to systems, pulse sequences,signal and diagnostic processors, diagnostic displays, and relatedmethods using novel application of nuclear magnetic resonance, includingmagnetic resonance spectroscopy, for diagnosing pain such as low backpain associated with degenerative disc disease.

Description of the Related Art

While significant effort has been directed toward improving treatmentsfor discogenic back pain, relatively little has been done to improve thediagnosis of painful discs—at least until quite recently.

Magnetic resonance imaging (MRI) is the primary standard of diagnosticcare for back pain. An estimated ten million MRIs are done each year forspine, which is the single largest category of all MRIs at an estimated26% of all MRIs performed. MRI in the context of back pain is sensitiveto changes in disc and endplate hydration and structural morphology, andoften yields clinically relevant diagnoses such as in setting ofspondylolysthesis and disc herniations with nerve root impingement (e.g.sciatica). In particular context of axial back pain, MRI is principallyuseful for indicating degree of disc degeneration. However, degree ofdisc degeneration has not been well correlated to pain. In one regard,people free of back pain often have disc degeneration profiles similarto those of people with chronic, severe axial back pain. In general, notall degenerative discs are painful, and not all painful discs aredegenerative. Accordingly, the structural information provided bystandard MRI exams of the lumbar spine is not generally useful fordifferentiating between painful and non-painful degenerative discs inthe region as related to chronic, severe back pain.

Accordingly, a second line diagnostic exam called “provocativediscography” (PD) is often performed after MRI exams in order tolocalize painful discs. This approach uses a needle injection ofpressurized dye in awake patients in order to intentionally provokepain. The patient's subjective reporting of pain level experiencedduring the injection, on increasing scale of 0-10, and concordancy tousual sensation of pain, is the primary diagnostic data used todetermine diagnosis as a “positive discogram”—indicating painfuldisc—versus a “negative discogram” for a disc indicating it is not asource of the patient's chronic, severe back pain. This has significantlimitations including invasiveness, pain, risks of disc damage,subjectivity, and lack of standardization of technique. PD has beenparticularly challenged for high “false+” rates alleged in variousstudies, although recent developments in the technique and studiesrelated thereto have alleged improved specificity of above 90%. (Wolferet al., Pain Physician 2008: 11:513-538, ISSN 1533-3159). However, thesignificant patient morbidity of the needle-based invasive procedure isnon-trivial, as the procedure itself causes severe pain and furthercompromises time from work. Furthermore, in another recent study PD wasshown to cause significant adverse effects to long term disc health,including significantly accelerating disc degeneration and herniationrates (on the lateral side of needle puncture). (Carragee et al., SPINEVolume 34, Number 21, pp. 2338-2345, 2009). Controversies around PDremain, and in many regards are only growing, despite the on-goingprevalence of the invasive, painful, subjective, harmful approach as thesecondary standard of care following MRI. PD is performed an estimated400,000 times annually world-wide, at an estimated total economic costthat exceeds $750 Million annually. The need for a non-invasive,painless, objective, non-significant risk, more efficient andcost-effective test to locate painful intervertebral discs of chronic,severe low back pain patients is urgent and growing.

A non-invasive radiographic technique to accurately differentiatebetween discs that are painful and non-painful may offer significantguidance in directing treatments and developing an evidence-basedapproach to the care of patients with lumbar degenerative disc disease(DDD).

Magnetic resonance spectroscopy (MRS), and related applications andpost-processing techniques, have been previously described. Thisincludes more recently in relation to intervertebral discs and certaindiagnostic applications related to medical conditions such asdegenerative disc disease and discogenic low back pain. However, suchprior efforts and disclosures, despite their alleged benefits, havenonetheless left open certain remaining needs and opportunities forvaluable new improvements still yet to be provided.

A need and opportunity still remain for improved MRS acquisition and/orpost-processing approaches for providing robust spectra representativeof chemical environments of tissues being analyzed. In particular, sucha need and opportunity exists to reduce risk of artifact contributionthat might confound chemical analysis of the tissue based upon processedMRS spectral data.

A need and opportunity also still remain for improved MRS approaches fordisc tissue analysis and diagnosis, such as for example for assisting inthe diagnosis of degenerative disc disease and/or discogenic low backpain.

SUMMARY OF THE INVENTION

One aspect of the present disclosure is an MRS post-processor providingenhanced frequency shift artifact correction in signal post-processingof acquired MRS data from tissue.

One aspect of the present disclosure is a MRS pulse sequence configuredto generate and acquire a diagnostically useful MRS spectrum from avoxel located principally within an intervertebral disc of a patient.

Another aspect of the present disclosure is an MRS signal processor thatis configured to select a sub-set of multiple channel acquisitionsreceived contemporaneously from multiple parallel acquisition channels,respectively, of a multi-channel detector assembly during arepetitive-frame MRS pulse sequence series conducted on a region ofinterest within a body of a subject.

Another aspect of the present disclosure is an MRS signal processorcomprising a phase shift corrector configured to recognize and correctphase shifting within a repetitive multi-frame acquisition seriesacquired by a multi-channel detector assembly during an MRS pulsesequence series conducted on a region of interest within a body of asubject.

Another aspect of the present disclosure is a MRS signal processorcomprising a frequency shift corrector configured to recognize andcorrect frequency shifting between multiple acquisition frames of arepetitive multi-frame acquisition series acquired within an acquisitiondetector channel of a multi-channel detector assembly during a MRS pulsesequence series conducted on a region of interest within a body of asubject.

According to one mode of this aspect, the frequency shift correctorcomprises a frequency shift estimator configured to estimate a frequencyshift error for the frames among the series, and a frequency shiftcorrector that corrects the frequency shift error based upon thefrequency shift estimation for the frames.

According to one embodiment of this mode, the frequency shift estimatoris configured to estimate the frequency shift error of the MRS data forthe frames in the time domain. According to another embodiment, thefrequency shift estimator is configured to estimate the frequency shifterror of the MRS data for the frames in the frequency domain. Accordingto another embodiment, the frequency shift corrector is configured tocorrect in the time domain the frequency shift error estimated by thefrequency shift estimator. According to another embodiment, thefrequency shift corrector is configured to correct in the frequencydomain the frequency shift error estimated by the frequency shiftestimator. According to another embodiment, both the frequency shiftestimator is configured to estimate the frequency shift error and thefrequency shift corrector is configured to correct the frequency shifterror in the time domain. According to another embodiment, both thefrequency shift estimator is configured to estimate the frequency shifterror and the frequency shift corrector is configured to correct thefrequency shift error in the frequency domain. According to anotherembodiment, the frequency shift estimator is configured to estimate thefrequency shift error in the time domain, and the frequency shiftcorrector is configured to correct the frequency shift error in thefrequency domain. According to another embodiment, the frequency shiftestimator is configured to estimate the frequency shift error in thefrequency domain, and the frequency shift corrector is configured tocorrect the frequency shift error in the time domain.

According to another embodiment, the frequency shift estimator isconfigured to estimate the frequency shift by comparing each frame witha reference by means of spectral cross-correlation. In one furtherembodiment, the reference is derived from multiple frames.

Another aspect of the present disclosure is a MRS signal processorcomprising a frame editor configured to recognize at least one poorquality acquisition frame, as determined against at least one thresholdcriterion, within an acquisition channel of a repetitive multi-frameacquisition series received from a multi-channel detector assemblyduring a MRS pulse sequence series conducted on a region of interestwithin a body of a subject.

According to one mode of this aspect, the frame editor is furtherconfigured to exclude said frames so recognized as poor relative qualityfrom the multi-frame acquisition series for further processing by one ormore additional signal processors. According to one embodiment of thismode, the further processing comprises frequency shift correction by afrequency shift corrector.

According to another mode of this aspect, the frame editor comprises awater peak confidence estimator configured to estimate a level or degreeof confidence in identifying a water peak signal in the spectral frame,and a water confidence threshold criterion. The frame editor is furtherconfigured to recognize and exclude frames from the series having anestimated degree of confidence for the water peak that do not meet thethreshold criteria. According to one embodiment, the degree ofconfidence is estimated as a percent (%) confidence level, on a scalefrom 0 to 100. According to a further embodiment, the thresholdcriterion comprises at least about 70 percent confidence. According toanother embodiment, the threshold criterion comprises at least about 90percent confidence. According to another embodiment, the degree ofconfidence is based at least in part upon an amplitude of the watersignal. According to another embodiment, the degree of confidence isbased at least in part upon a full width half max amplitude (FWHM)measurement of the water signal.

According to another mode of this aspect, the frame editor comprises afull width half maximum (FWHM) test module, which comprises a FWHMmeasurement of the water peak and a FWHM measurement thresholdcriterion. The frame editor is further configured to recognize andexclude frames from the series having a computed FWHM of the water peakthat exceeds the threshold criteria.

According to another mode of this aspect, the frame editor comprises afrequency error test module, comprising a frequency error measurement ofthe water peak and a frequency error threshold criterion. The frequencyerror is the difference between the computed location of the water peak(maximum peak of the spectrum) and the defined location of the waterpeak (DC). The editor is further configured to recognize and excludeframes from the series having a frequency error that exceeds thethreshold criteria.

Another aspect of the present disclosure is an MRS signal processor thatcomprises an apodizer to reduce the truncation effect on the sampledata. The apodizer can be configured to apodize an MRS acquisition framein the time domain otherwise generated and acquired by via an MRS aspectotherwise herein disclosed, and/or signal processed by one or more ofthe various MRS signal processor aspects also otherwise hereindisclosed.

Another aspect of the present disclosure is an MRS diagnostic processorconfigured to process information extracted from an MRS spectrum for aregion of interest in a body of a subject, and to provide the processedinformation in a manner that is useful for diagnosing a medicalcondition or chemical environment associated with the region ofinterest.

Another aspect of the present disclosure is an MRS system comprising anMRS pulse sequence, MRS signal processor, and MRS diagnostic processor,and which is configured to generate, acquire, and process an MRSspectrum representative of a region of interest in a body of a patientfor providing diagnostically useful information associated with theregion of interest.

Another aspect of the present disclosure is an MRS system which adjustsa spectrum in an upfield range relative to water (e.g. <about 4.7 ppm)based upon a measured parameter in a downfield range relative to water(e.g. >about 4.7 ppm). According to one mode, at least one region of anupfield range of less than about 3.5 ppm is adjusted based upon at leastone roughly mirrored opposite range of greater than about 5.9 or 6 ppm.These may be accomplished within about +/−0.5 ppm. According to oneembodiment, the downfield range is reflected around a central line tomirror the upfield range, and is subtracted from the upfield range.According to one further embodiment, the central line is at about water(e.g. about 4.7 ppm). According to another embodiment, the downfieldrange is further aligned following reflection for enhanced correlationbetween the downfield range signals and the upfield range signals.According to another embodiment, the correlation between the aligneddownfield and upfield regions is measured and used as a test to apply orbypass the adjustment. According to another embodiment, a signal:noiseratio (SNR) calculation of the adjusted upfield spectrum followingsubtraction is compared with an SNR calculation prior to the adjustmentin order to determine whether to use the adjusted spectrum or discard itand retain the unadjusted spectrum According to another embodiment, thecorrelation between the reflected and aligned downfield spectrum and theupfield spectrum ranges is used to determine and identify a degree ofpotential artifact in the upfield spectrum signals.

Still further aspects of the present disclosure comprise various MRSsystem and method aspects associated with the other MRS system,sequence, and processor aspects described above.

For example, another such aspect comprises processing MRS data from anacquisition series conducted across first and second phase groupscomprising first and second sets of acquisition frames acquired at firstand second different respective phases along a phase cycle. Further tothis aspect, the system and method processes at least in part the firstand second phase cycle groups separately. According to one mode of thisaspect, the separate phase group processing comprises at least one ofthe other processing aspects, modes, embodiments, variations, andfeatures elsewhere herein described, such as for example one or more offrequency error estimation, frame editing, frequency correction, andphase correction. According to another mode, the first and second set ofacquisition frames for each phase group are combined into one phasegroup combined result. According to another mode, the first and secondphase groups are combined to provide an averaged spectral resultfollowing the separate processing. In another mode, the first and secondsets of each phase group are combined, and then the combined phasegroups are combined. Another mode comprises reducing or removingartifact separately between the phase groups. In another mode, theartifact is removed by combining the phase groups after separateprocessing and in-group frame combining. In another mode, frame editingis performed on each phase group, and a different number of frames areretained (or conversely filtered out) in each group prior to combiningin-group frames and the groups together.

Another aspect of the present disclosure is a method for processing aset of multiple serially acquired magnetic resonance spectroscopy (MRS)free induction decay (FID) frames from a multi-frame MRS acquisitionseries from a region of interest (ROI) in a subject, and for providing apost-processed MRS spectrum. This method comprises:

dynamically varying a parameter of a processing step between multipleparameter values affecting a corresponding change between multiplefeature values of a feature of a processed result of performing the stepon the set;

performing the step on the set at each of the multiple parameter valuesand providing multiple processed results, respectively, comprisingmultiple respective feature values, respectively;

comparing the feature values of the processed results;

determining a chosen feature value among the multiple feature valuesbased upon the comparison and corresponding to a chosen parameter value;and

setting the parameter to the chosen parameter value for performing theprocessing step to provide the post-processed MRS spectrum.

According to one mode of this method, the chosen feature value comprisesa best feature value relative to a quality criterion for the value.

According to another mode, the feature comprises signal:noise ratio(SNR) of a spectral peak region of the post-processed MRS spectrum.

According to another mode, the feature comprises a line width of aspectral peak region of the post-processed MRS spectrum.

According to one embodiment of this mode, the line width comprises afull width half max (FWHM) measurement of a spectral peak region of thepost-processed MRS spectrum.

According to another mode, the feature comprises signal:noise ratio(SNR) or full width half max (FWHM) of a spectral peak region of thepost-processed MRS spectrum; and the spectral peak region comprises atleast one of a water region and a second region that is different thanthe water region and corresponding with a chemical of interest for theMRS acquisition.

According to one embodiment of this mode, the second region correspondswith an n-acetyl acetate (NAA) chemical bond.

According to another mode, the processing step comprises: frame editingby comparing a quality value for a quality of each frame against athreshold value for the quality and determining if the frame isqualified and included, or unqualified and excluded, in the set forfurther processing; and the parameter comprises the threshold value,such the multiple parameter values comprise multiple threshold values,and such that varying the multiple threshold values corresponds withvarying the frames which are qualified and included, or unqualified andexcluded, from the set for further processing in providing thepost-processed MRS spectrum.

According to one further embodiment, the quality comprises a confidenceestimating a water peak location.

According to another further embodiment, the quality comprises afrequency error of a water peak location relative to a referencefrequency.

According to another further embodiment, the quality comprises anamplitude of a water peak.

According to another further embodiment, the quality comprises a fullwidth half max (FWHM) of a water peak.

According to another mode, the processing step comprises frame editingby comparing a quality value for a quality of each frame against athreshold value for the quality and determining if the frame isqualified and included, or unqualified and excluded, in the set forfurther processing, and also comprises bypassing frame editing if anumber of qualified frames is below a threshold number; and theparameter comprises the threshold number.

Another mode comprises: dynamically varying multiple said parameters ofthe processing step between multiple said parameter values according toa multi-variate matrix comprising multiple combinations of the variedparameter values for each said parameter; performing the step on the setat each combination of parameter values and providing said multipleprocessed results, respectively, comprising said multiple respectivefeature values, respectively; determining the chosen feature value amongthe multiple feature values based upon the comparison and correspondingto a chosen combination; and providing the post-processed MRS spectrumby running the step on the set with the multiple parameters assigned tothe chosen combination of respective parameter values.

According to one embodiment of this mode, the multiple parameterscomprise multiple frame editing criteria for qualifying and including,or unqualifying and excluding, frames from the set for furtherprocessing steps to provide the post-processed MRS spectrum.

Another aspect of the present disclosure is also a method for processinga set of multiple serially acquired magnetic resonance spectroscopy(MRS) free induction decay (FID) frames from a multi-frame MRSacquisition series from a region of interest (ROI) in a subject, and forproviding a post-processed MRS spectrum. The method according to thisaspect comprises:

providing a first post-processed MRS spectrum;

identifying a feature of a downfield region of the spectrum that islocated downfield along a larger part per million (ppm) range than awater peak location;

comparing the feature with a reference;

performing a processing operation based on the comparison and related tothe post-processed MRS spectrum.

According to one mode of this aspect, the identified feature comprises adownfield peak region of the downfield region.

According to another mode, the comparison comprises comparing thedownfield peak region against an upfield peak region of an upfieldregion of the spectrum located along a smaller ppm range than a waterpeak location.

According to another mode, the processing operation comprises adjustingthe upfield region based on the comparison.

According to one embodiment of this mode, the adjusting comprisesreducing at least a portion of the upfield region by a subtractionoperation based upon the comparison.

According to another mode, the comparison comprises mapping a reflectionof the downfield region around a center line against the upfield region.

According to one embodiment of this mode, the center line is located atabout the water peak location.

According to another embodiment, the comparison further comprisesadjusting the center line to an offset relative to the water peaklocation to correspond with an optimum correlation between at least aportion of the reflected downfield region and a corresponding portion ofthe upfield region:

According to a further embodiment, the center line adjustment isperformed within a range of offset.

Another embodiment comprises baseline correcting the reflected downfieldand upfield regions, respectively, and subtracting the baselinecorrected reflected downfield region from the baseline corrected upfieldregion.

According to another mode of the present aspect, the processingoperation comprises qualifying the post-processed MRS spectrum basedupon the identified feature.

Another mode comprises adjusting a water suppression parameter of theMRS acquisition based upon the comparison.

Another aspect of the present disclosure is also a method for processinga set of multiple serially acquired magnetic resonance spectroscopy(MRS) free induction decay (FID) frames from a multi-frame MRSacquisition series from a region of interest (ROI) in a subject, and forproviding a post-processed MRS spectrum. The method according to thisaspect comprises:

determining a frequency shift error of an FID frame by performingspectral cross-correlation between the FID frame and a reference frame.

One mode of this aspect comprises shifting the FID frame by thefrequency shift error so as to increase coherence of the FID frame withthe reference spectrum.

Another mode comprises forming the reference spectrum by averagingmultiple FID frames from the acquisition series.

Another mode comprises: performing an initial frequency shift errorcorrection operation by estimating a frequency shift error for each FIDframe of the series by determining a location of a water peak bylocating maximum peak value in a range around an expected location forthe water peak, calculating a difference between the determined locationand the expected location, and adjusting the FID frame by thedifference; and averaging the respectively shifted FID frames to formthe reference spectrum.

According to another mode, the determining is performed on each FIDframe of the series.

According to another mode, the series comprises first and second groupsof FID frames, the FID frames of the first group comprise a first phase,and the FID frames of the second group comprises a second phase, and thedetermining is performed on each FID frame for each group separatelyfrom the other group.

According to another mode, the reference frame comprises another FIDframe from the series.

Another mode comprises: comparing the frequency shift error against athreshold value; and qualifying and including, or disqualifying andexcluding, the FID frame from the set for further processing to providethe post-processed MRS spectrum based upon the comparison.

Another aspect of the present disclosure is also a method for processinga set of multiple serially acquired magnetic resonance spectroscopy(MRS) free induction decay (FID) frames from a multi-frame MRSacquisition series from a region of interest (ROI) in a subject, and forproviding a post-processed MRS spectrum. The method according to thisaspect comprises adjusting a phase of an FID frame by a phase shiftconfigured to reduce an extent of the absorption spectrum having anegative value.

One mode of this aspect further comprises performing at least one of agolden section search and a parabolic interpretation to compute a fastfourier transform (FFT) on the FID frame within an iteration loop.

Another aspect of the present disclosure is also a method for processinga set of multiple serially acquired magnetic resonance spectroscopy(MRS) free induction decay (FID) frames from a multi-frame MRSacquisition series from a region of interest (ROI) in a subject, and forproviding a post-processed MRS spectrum. The method according to thisaspect comprises: estimating a line width between two opposite walls ofa water peak region of an FID frame of the series; comparing theestimated line width against a threshold value; and qualifying andincluding, or disqualifying and excluding, the FID frame from the setfor further processing to provide the post-processed MRS spectrum basedat least in part upon the comparison.

According to one mode of this aspect, the line width comprises a fullwidth half max (FWHM) of the water peak region.

Another mode comprises: calculating an average FWHM for the included FIDframes; comparing the average FWHM against an average threshold value;and assigning a quality indicator to the post-processed MRS spectrumbased upon the comparison.

Another mode comprises adjusting a water suppression parameter basedupon the comparison.

According to another mode, the water suppression parameter compriseswater suppression bandwidth.

According to another mode, the water suppression parameter comprises adegree of water suppression.

Another mode further comprises determining the threshold value basedupon a tesla strength of an MRS system from which the MRS acquisitionseries was acquired.

Another aspect of the current disclosure is also a method for processinga set of multiple serially acquired magnetic resonance spectroscopy(MRS) free induction decay (FID) frames from a multi-frame MRSacquisition series from a region of interest (ROI) in a subject, and forproviding a post-processed MRS spectrum. The method according to thisaspect comprises:

associating first and second groups of FID frames within the set withfirst and second different relative phases, respectively, of a phasestep cycle along the MRS acquisition series;

processing the first and second groups separately;

averaging the first and second groups separately to form first andsecond average interim spectra; and

averaging the first and second average interim spectra.

One mode of this aspect comprises performing at least one of phase shiftcorrection, frame editing, and frequency shift correction separately forthe first and second groups.

Another mode comprises performing the method on a number equal to Ngroups of FID frames, the FID frames comprising similar phasing withineach group but different than in other groups, processing each of the Ngroups separately, averaging the FID frames separately within eachgroup, averaging the group averages together, and wherein N equals aninteger more than two.

According to one embodiment of this mode, the MRS acquisition seriescomprises a number F of FID frames, the different relative phasescomprise a unique number S of phase steps, and F divided by S equals N,wherein N comprises a whole number.

Another aspect of the present disclosure is also a method for processinga set of multiple serially acquired magnetic resonance spectroscopy(MRS) free induction decay (FID) frames from a multi-frame MRSacquisition series from a region of interest (ROI) in a subject, and forproviding a post-processed MRS spectrum. The method according to thisaspect comprises:

generating a first baseline estimate using a rank order filter (ROF) ona range of interest of an interim MRS spectrum partially post-processedfrom the set;

fitting a polynomial to at least a portion of the first baselineestimate to generate a baseline curve; and

subtracting the baseline curve from the interim MRS spectrum to providethe post-processed MRS spectrum.

One mode of this aspect also comprises determining first and secondregions of the first baseline estimate to be included and excluded,respectively, for the polynomial fit; and fitting the polynomial to onlythe first regions of the baseline estimate.

Another aspect of the present disclosure is a method performed accordingto any one of claims 1 to 50 using a processor.

One mode of this aspect further comprises controlling the processor toperform said one or more of said methods using a computer program storedin a non-transitory computer readable medium.

Another aspect of the present disclosure comprises performing one ormore of the methods elsewhere herein described, wherein: the ROIcomprises a portion of a musculoskeletal joint.

According to one mode of this aspect, the joint comprises a skeletaljoint.

According to one embodiment of this mode, the skeletal joint comprises aspinal joint.

According to another further embodiment, the ROI comprises at least aportion of an intervertebral disc.

Another aspect of this disclosure comprises performing any one or moreof the methods elsewhere described herein, wherein said MRS acquisitionseries is acquired from a single voxel positioned within the ROIaccording to a single voxel spectroscopy pulse sequence

Another aspect of this disclosure also comprises performing any one ormore of the methods elsewhere described herein, wherein said MRSacquisition comprises a multi-voxel spectroscopy acquisition comprisingmultiple said MRS acquisition series corresponding to each said voxel,and performing the method on each said MRS acquisition series.

Another aspect of the present disclosure also comprises performing anyone or more of the methods described elsewhere herein, and furthercomprising: providing an MRS scanner system; and acquiring the MRSacquisition series using the MRS scanner system on the subject.

Another aspect of the present disclosure comprises one or morenon-transitory computer readable media comprising computer instructionsconfigured to cause one or more computer processors to perform actionscomprising one or more of the methods described elsewhere herein:

Another aspect of the present disclosure comprises one or morenon-transitory computer readable media comprising computer instructionsconfigured to cause one or more computer processors to perform actionscomprising one or more of the methods described elsewhere herein.

Another aspect of the present disclosure is a method for processing aset of multiple serially acquired magnetic resonance spectroscopy (MRS)free induction decay (FID) frames from a multi-frame MRS acquisitionseries from a region of interest (ROI) in a subject, and for providing apost-processed MRS spectrum. The method according to this aspectcomprises performing a combination of at least two of the following:

performing the step on the set at each of the multiple parameter valuesand providing multiple processed results, respectively, comprisingmultiple respective feature values, respectively, comparing the featurevalues of the processed results, determining a chosen feature valueamong the multiple feature values based upon the comparison andcorresponding to a chosen parameter value, setting the parameter to thechosen parameter value for performing the processing step to provide thepost-processed MRS spectrum;

providing a first post-processed MRS spectrum, identifying a feature ofa downfield region of the spectrum that is located downfield along alarger part per million (ppm) range than a water peak location,comparing the feature with a reference, performing a processingoperation based on the comparison and related to the post-processed MRSspectrum;

determining a frequency shift error of an FID frame by performingspectral cross-correlation between the FID frame and a reference frame;

adjusting a phase of an FID frame by a phase shift configured to reducean extent of the absorption spectrum having a negative value.

estimating a line width between two opposite walls of a water peakregion of an FID frame of the series, comparing the estimated line widthagainst a threshold value, and qualifying and including, ordisqualifying and excluding, the FID frame from the set for furtherprocessing to provide the post-processed MRS spectrum based at least inpart upon the comparison;

associating first and second groups of FID frames within the set withfirst and second different relative phases, respectively, of a phasestep cycle along the MRS acquisition series, processing the first andsecond groups separately, averaging the first and second groupsseparately to form first and second average interim spectra, andaveraging the first and second average interim spectra; and

generating a first baseline estimate using a rank order filter (ROF) ona range of interest of an interim MRS spectrum partially post-processedfrom the set, fitting a polynomial to at least a portion of the firstbaseline estimate to generate a baseline curve, and subtracting thebaseline curve from the interim MRS spectrum to provide thepost-processed MRS spectrum.

Another aspect of the present disclosure comprises one or morenon-transitory computer readable media comprising computer instructionsconfigured to cause one or more computer processors to perform actionscomprising one or more of the methods described elsewhere herein:

One further mode of each of the various aspects noted above, andcombinations therebetween, comprises a magnetic resonance (MR) systemcomprising a magnet. According to one embodiment, a controller isprovided for controlling the magnet to conduct an MR exam of a patient.According to another further embodiment, the controller is configured torun a pulse sequence to actuate the magnetic to induce a pulsatilemagnetic field in an area corresponding with a patient or tissue, suchas according to certain pulse sequence aspects disclosed herein and/orother pulse sequence embodiments as may be contemplated by one ofordinary skill to the extent consistent with this disclosure. Further tothis embodiment, the pulsatile magnetic field induced by the pulsesequence is configured to invoke a tissue response, unique to theinduced pulsatile magnetic field, in a region of interest of the tissuedesignated for evaluation. The tissue response uniquely invoked by theinduced pulsatile magnetic field generates signals along a frequencyspectrum which emanate from the region of interest and are captured byan antenna receiver coil also provided in spatial relation to the regionof interest. The signals of the invoked tissue response comprise uniquefrequency components corresponding with unique respective chemicalconstituents in the region of interest of the tissue. Accordingly, thespectrum of frequency response components invoked in the region ofinterest provides MRS information related to the tissue chemistry inthat region. In a further embodiment, the region of interest comprisesat least a portion of an intervertebral disc. In still a furtherembodiment, a voxel is prescribed to correspond with the region ofinterest and used to determine where spatially the induced and acquiredMRS spectrum information is to correspond.

According to another further mode of each of the various aspects notedabove, a non-transitory computer readable medium encoded with a computerprogram, and configured to be run by a processor, is provided. Accordingto one embodiment of this mode, a processor configured to run thecomputer program is also provided. According to still a furtherembodiment, the program run by the processor comprises at least one ofthe MRS signal processing aspects herein disclosed, such as for examplechannel selection, phase correction, frame editing, frequency shiftcorrection, apodization, baseline correction, downstream artifactcorrection, metabolite range measurement, and diagnostic interpretativeprocessing.

Each of the foregoing aspects, modes, embodiments, variations, andfeatures noted above, and those noted elsewhere herein, is considered torepresent independent value for beneficial use, including even if onlyfor the purpose of providing further combination with others, andwhereas their various combinations and sub-combinations are furthercontemplated aspects also of independent value for beneficial use, asmay be made by one of ordinary skill based upon a thorough review ofthis disclosure in its entirety.

Other aspects not specifically described above are also contemplated asmade clear in the detailed description below. For example, additionalaspects of the present disclosure comprises respective methodscorresponding with manufacturing and using the systems and devicesdescribed in these aspects above, and in the description below—as wouldbe apparent to one of ordinary skill. In addition, further more detailedmodes, embodiments, features, and variations of such aspects describedabove and below are also herein contemplated.

BRIEF DESCRIPTION OF THE DRAWINGS

These and other features, aspects, and advantages of the presentdisclosure will now be described with reference to the drawings ofembodiments, which embodiments are intended to illustrate and not tolimit the disclosure.

FIGS. 1-43C are variously described in the detailed description providedbelow.

FIG. 1 shows a schematic flow chart for certain steps of a method aspectof the current disclosure.

FIG. 2 shows a schematic flow chart for certain steps of another methodaspect of the disclosure.

FIG. 3 shows a schematic flow chart for certain steps of another methodaspect of the disclosure.

FIG. 4 shows a diagrammatic illustration of a peak and with respect tocertain features thereof, as reference to certain aspects of the presentdisclosure.

FIG. 5 shows a residual water full bandwidth (BW) half amplitude byframe, also referred to as full width half max (FWHM), for a 192 frameacquisition, and also shows a frame editing threshold such that frameswith FWHM above the threshold are excluded.

FIG. 6A shows another residual water full bandwidth half amplitude byframe plot over another 192 frame MRS acquisition from a lumbar disc,and shows a change around 80 frames into the acquisition where someframes begin demonstrating elevated FWHM bandwidth relative to a secondhigher threshold criteria versus that used reflected in the FIG. 5example.

FIG. 6B shows a time-intensity plot vs. frequency for each frame in a192 frame acquisition prior to frequency correction, where the y-axis isframes in acquisition ordered numbering (e.g. 1 to 192) moving downwardvertically, and frequency range horizontally on y-axis, with eachacquired spectral frame reflected as single horizontal line withspectral amplitudes reflected by contrast (brightness) value setting(higher amplitudes=lighter contrast).

FIG. 6C shows the same time-intensity plot for the same serial frame MRSacquisition shown in FIG. 6A, except after frequency shift correction onretained qualified frames performed according to one embodiment.

FIG. 6D shows the spectrum for the same acquisition featured in FIGS.6A-C, after phase and frequency adjustment, according to one embodimentof the disclosure.

FIG. 7A shows a residual water full BW half amplitude by frame plot forthe same serial frame MRS acquisition in the lumbar disc of FIGS. 6A-D,similar to the plot shown in FIG. 6A, but in a process according toanother embodiment applying a different lower threshold criteria to editout more, and retain less qualified, frames.

FIG. 7B shows a similar time-intensity plot to that shown in FIG. 6B,but for the acquisition shown in FIG. 7A.

FIG. 7C shows a similar time-intensity plot to that shown in FIG. 6C,but for the acquisition shown in FIG. 7A and after frame editing andfrequency correction.

FIG. 7D shows the spectrum for the same acquisition featured in FIGS.7A-C, after frame editing and phase and frequency adjustment of retainedqualified frames, according to another example embodiment.

FIG. 8A shows a plot for confidence of frequency estimate for each framealong the same 192 frame serial MRS acquisition featured in FIGS. 6A-7D,against which a threshold value of 0.5 (not shown) is applied such thatframes must be below the threshold to be qualified and retained forfurther processing (unless editing is bypassed).

FIG. 8B shows a similar time-intensity plot to that shown in FIG. 7B,but for the acquisition shown in FIG. 8A.

FIG. 8C shows a similar time-intensity plot to that shown in FIG. 7C,but for the acquisition shown in FIG. 8A and after frame editing andfrequency correction according to the thresholds and operatingparameters of the example featured in this series of FIGS. 8A-D.

FIG. 8D shows the spectrum for the same acquisition featured in FIG. 8C,after frame editing and phase and frequency adjustment of retainedqualified frames, according to another example embodiment.

FIG. 9A shows a plot for confidence of frequency estimate for each framesimilar to that shown in FIG. 8A and for the same acquisition featuredin FIGS. 6A-8D, but with a threshold criteria of 0.2 (not shown).

FIG. 9B shows a similar time-intensity plot to that shown in FIG. 8B.

FIG. 9C shows a similar time-intensity plot to that shown in FIG. 9C,and after frame editing and frequency correction according to thethresholds and operating parameters of the example featured in thisseries of FIGS. 9A-D.

FIG. 9D shows the spectrum similar to that shown in FIG. 8D, butreflecting the retained frames shown in FIG. 9C after frame editing andphase and frequency adjustment of retained qualified frames according tothe different threshold criteria of this example.

FIG. 10 shows a post-processed spectrum from a lumbar disc spectroscopyexam, prior to baseline correction, in full spectrum range view.

FIG. 11 shows an input spectrum according to an illustrated downfieldartifact correction embodiment of Example 1, and shows water peak andsuppression side-lobe artifacts marked.

FIG. 12A shows the same input spectrum as FIG. 11, but with addition ofa mirrored spectrum in overlay, full range including up and downfieldspectral regions.

FIG. 12B shows the same input spectrum and mirrored spectrum of FIG.12A, with offset downfield correction by 1 bin, in zoomed metaboliterange.

FIG. 13A shows the same spectrum of Example 1 shown in FIGS. 10-12B, butwith the spectrum shifted to water and mirrored spectrum, full range.

FIG. 13B shows again the same spectrum of Example 1 shown in FIG. 13A,zoomed to metabolite range of interest below 4.5 ppm on chemical shiftspectrum.

FIG. 14A shows the spectral processing in further steps under Example 1,but with fine alignment having been performed between the spectrum andmirrored spectrum, shown in full range.

FIG. 14B shows same interim result of FIG. 14A, but in zoomed view alongthe metabolite range of interest.

FIG. 15 shows a further progressed stage of Example 1, showing thespectrum and baseline corrected spectrum in overlay, full spectrumrange.

FIG. 16 shows the mirrored spectrum in overlay with baseline correctedmirrored spectrum of Example 1, full spectrum range.

FIG. 17 shows aligned baseline corrected up-field and mirror down-fieldspectra, in zoomed metabolite range.

FIG. 18 shows adjusted up-field spectrum using baseline correctedspectrum, zoomed metabolite range.

FIG. 19 shows aligned up-field and mirror down-field spectra in overlayaccording to a further embodiment under Example 1, without baselinecorrection, in zoomed metabolite range.

FIG. 20A shows an adjusted up-field spectrum using non-baselinecorrected spectra of FIG. 19A, zoomed metabolite range.

FIG. 20B shows an adjusted spectrum post-fine alignment and subtraction,baseline corrected, per a further embodiment under Example 1, zoomedmetabolite range.

FIG. 21A shows an input spectrum according to another downfieldcorrection embodiment of Example 2, post-processed (but without baselinecorrection or downfield adjustment), in full spectrum range.

FIG. 21B shows the input spectrum of FIG. 21A further to the Example 2embodiment, in zoomed metabolite spectrum range.

FIG. 22A shows another graph illustrating the Example 2 embodiment shownin FIGS. 21A-B, but shows a reflected or mirrored representation of thedownfield spectrum (mirrored around water) in overlay over thecorresponding opposite upfield range of the same spectrum, followingcourse alignment, in full spectrum range.

FIG. 22B shows the same spectrum and reflected overlay downfieldspectrum shown in FIG. 22A further to Example 2, in zoomed metaboliterange.

FIG. 23A shows the same overlay view of upfield spectrum with reflecteddownfield spectrum of Example 2, but post fine alignment, in fullspectrum range.

FIG. 23B shows the same spectral overlay of Example 2 shown in FIG. 23A,but in zoomed metabolite spectrum.

FIG. 24A shows the same overlay view of spectra under Example 2 shown inFIG. 23B, but after baseline correction, in zoomed spectrum range.

FIG. 24B shows the result under Example 2 of subtracting and correctingthe spectrum of FIG. 24A by subtraction with the reflected downfieldspectrum also shown in FIG. 24A, baseline corrected, in zoomedmetabolite range.

FIG. 25A shows another graph under Example 2 showing overlay betweenreflected downfield spectrum onto upfield spectrum, post-fine alignmentand without baseline correction, in zoomed spectrum range.

FIG. 25B shows adjusted spectrum after subtracting reflected baselinefrom upfield spectra shown in FIG. 25A, without baseline correction, inzoomed metabolite spectrum range.

FIG. 25C shows the adjusted spectrum of FIG. 25B, post-baselinecorrection, in zoomed metabolite spectrum range.

FIG. 26A shows an input spectrum according to another embodiment ofExample 3, and which is post-processed (but without baseline correctionor downfield adjustment), in full spectrum range.

FIG. 26B shows the input spectrum of Example 3 shown in FIG. 26A, inzoomed spectrum range.

FIG. 27A shows the same spectrum of Example 3 shown in FIG. 26A, but inoverlay with a reflected spectrum mirrored around water signal, withcoarse alignment, in full spectrum range.

FIG. 27B shows the same overlaid spectra of Example 3 shown in FIG. 27A,in zoomed metabolite range.

FIG. 28A shows the same overlaid spectra of Example 3 shown in FIG. 27A,but after fine alignment adjustment along the chemical shift X-axis, infull spectrum range.

FIG. 28B shows the same overlaid spectra post-fine alignment of Example3 shown in FIG. 28A, in zoomed metabolite range.

FIG. 29A shows the overlaid spectra post-fine alignment under Example 3shown in FIG. 28B, after baseline correction, in zoomed metabolitespectrum range.

FIG. 29B shows an adjusted spectrum result under Example 3 aftersubtracting the fine aligned baseline corrected spectra shown in overlayin FIG. 29A, and after baseline correction, in zoomed metabolitespectrum range.

FIG. 30A shows overlay plot of the Example 3 spectrum and reflecteddownfield spectrum, without baseline correction, in zoomed metaboliterange.

FIG. 30B shows an adjusted spectrum derived after subtracting thereflected from the upfield spectra shown in FIG. 30A, post-finealignment but without baseline correction, in zoomed metabolite spectrumrange.

FIG. 30C shows the adjusted spectrum shown in FIG. 30B, but afterbaseline correction, in zoomed metabolite spectrum range.

FIG. 31A shows an input spectrum and reflected spectrum thereof mirroredabout water according to another Example 4, post processed (but withoutbaseline correction or downfield adjustment), full spectrum range.

FIG. 31B shows the same input spectrum and reflected spectrum in overlayas shown in FIG. 31A, but in zoomed metabolite range.

FIG. 32A shows the same overlaid spectra as shown in FIG. 31A, but aftercoarse alignment, in full spectrum range.

FIG. 32B shows the same spectral overlay after coarse alignment ofExample 4 shown in FIG. 32A, in zoomed metabolite spectrum range.

FIG. 33A shows the same overlaid spectra of Example 4 shown in FIG. 32A,after fine alignment, in full spectrum range.

FIG. 33B shows the same overlaid spectra of Example 4 shown in FIG. 33A,in zoomed metabolite range.

FIG. 34A shows the same overlaid spectra of Example 4 shown in FIG. 33A,each after baseline correction, in zoomed metabolite spectrum range.

FIG. 34B shows an adjusted spectrum of Example 4 after subtracting thefine-aligned, baseline spectra overlaid in FIG. 34A, zoomed metaboliterange.

FIG. 35A shows a graph overlay of reflected downfield vs. upfieldspectra of Example 4, post-fine alignment, but without baselinecorrection, in zoomed spectrum range.

FIG. 35B shows adjusted spectral results after subtracting the overlaidspectra shown in FIG. 35A of Example 4, without baseline correction,zoomed metabolite spectrum range.

FIG. 35C shows the same adjusted spectrum as shown in FIG. 35B under theExample 4, but after post-adjustment baseline correction, in zoomedmetabolite spectrum range.

FIG. 36 shows an x-y plot of input versus reflected input (“tupnI”) datapoints from an acquisition to validate the mirroring function (asmirrored around point or “bin” 9).

FIG. 37 shows another x-y plot of input versus reflected input datapoints from an acquisition for similar validation as demonstrated in theFIG. 36 plot, although showing a calculated reflection center equal toabout 12 in the particular example shown.

FIG. 38 shows a spectrum from a lumbar disc acquisition in overlaybetween two different aggregated phase groups from a 2-step phase cycledMRS acquisition series, and average between them.

FIG. 39 shows a flow diagram of an automated method to perform frequencyshift estimation and correction, and frame selection in the context ofphase group processing.

FIG. 40A shows the full spectrum of a single frame from a lumbar MRSdisc acquisition.

FIG. 40B shows the metabolite range for this same frame shown in

FIG. 40A.

FIG. 40C shows a spectrum representing the average of all 192 framescollected in the series featured in FIGS. 40A-B, and with only veryslight phase correction and apodization processing.

FIG. 40D shows this same spectrum shown in FIG. 40C, but in a zoomedview of a metabolite range of interest.

FIG. 41A shows a post-processed spectrum from an MRS acquisition from alumbar disc, baseline corrected using polynomial fit to an ROFincorporating all spectral data in the range (including outliers).

FIG. 41B shows certain baseline correction details overlaid to thespectrum, full spectrum range, for the acquisition featured in FIG. 41A,with the ROF is shown below the spectrum in the metabolite range, thebold dotted line showing polynomial ft to the ROF without outliers,versus the narrower dashed line that shows the polynomial fit withoutliers.

FIG. 41C shows these baseline correction details shown in FIG. 41B, in azoomed view more narrowly drawn to a metabolite range of interest fromthe spectrum.

FIG. 41D shows a final spectra post baseline correction using thepolyfit without outliers contributing to the ROF used for the fit.

FIG. 42A shows a full range spectrum from an intravertebral disc MRSacquisition, post-processed according to a combination of multipleembodiments of the current disclosure.

FIG. 42B shows a time-intensity plot of the individual frames along theacquisition series providing the post-processed spectrum shown in FIG.42A, and shows prior to processing (top), with frame editing (middle),and after frequency correction of the retained qualified frames (bottom)(unqualified screened out frames replaced with dark horizontal lines inthe plot).

FIG. 42C shows a frame editing panel comprising a series of water signalanalyses comprising per-frame plots of “frequency estimate confidence”(confidence of peak location), frequency error, in-phase amplitude,FWHM, and finally qualified retained vs. excluded unqualified (red)frames at bottom.

FIG. 43A shows a similar spectral plot for the same acquisition shown inFIG. 42A, but bypassing frequency correction by setting the FWHMthreshold criteria below the FWHM of all frames (e.g. excluding them allthus bypasses frame editing and frequency correction all together).

FIG. 43B shows a time intensity plot for the same series of phases ofthe same acquisition shown in FIG. 43B.

FIG. 43C shows a frame editing panel comprising the water signalanalyses similar to that shown in FIG. 42C, and on the same acquisitionfeatured in FIG. 42C, but related to the different frame editing andfrequency correction thresholds and related processing featured in FIGS.43A-B.

DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS

This disclosure relates to the following co-pending US patentapplications, which were previously filed and incorporated herein intheir entirety by reference thereto: US 2008/0039710, US 2009/00300308,US 2011/0087087 (U.S. Ser. No. 12/579,371 to Peacock et al., filed Oct.14, 2009).

This disclosure relates to the following co-pending Published PCT PatentApplications, which were previously filed and incorporated herein intheir entirety by reference thereto: WO 2006/081471, WO 2007/035906, WO2011/047197 (International Patent Application No. PCT/US2010/052737 toPeacock et al., filed 14 Oct. 2010).

Many alternative embodiments of the present aspects may be appropriateand are contemplated, including as described in these detailedembodiments, though also including alternatives that may not beexpressly shown or described herein but as obvious variants or obviouslycontemplated according to one of ordinary skill based on reviewing thetotality of this disclosure in combination with other availableinformation. For example, it is contemplated that features shown anddescribed with respect to one or more particular embodiments, may alsobe included in combination with another embodiment even though notexpressly shown and described in that specific combination.

For purpose of efficiency, where reference numbers may be used in theFigures, they may be repeated between the Figures where they areintended to represent similar features between otherwise variedembodiments, though those features may also incorporate certaindifferences between embodiments if and to the extent specified as suchor otherwise apparent to one of ordinary skill.

Various aspects of the present disclosure principally relates to amagnetic resonance spectroscopy (MRS) pulse sequence and/or relatedpost-processing system and method. One aspect accordingly combines thedetailed embodiments, either individually or in combination with one ormore others, into an overall system. One such system and methodcontemplated according to one aspect of this disclosure is shown inFIG. 1. This approach is considered highly beneficial and includeschannel selection via a channel selector, phase correction via a phasecorrector, frame editing via a frame editor, apodization via anapodizer, transformation via a transformer, frequency correction via afrequency corrector, and combined averaging via a combiner. In addition,further embodiments also include in one regard a baseline correction viaa baseline corrector, and in another regard reflected downfield artifactcorrection via another appropriately deployed corrector. An additionaldetailed embodiment is shown in FIG. 2, which further illustrates theflow of MRS signal processing according to these system components, andfurther downstream aspects applying post-processed spectral results tospectral measurement and diagnostic interpretive processing, asconsistent with that particular embodiment. Additional particularlybeneficial present embodiments of these various aspects of the presentdisclosure are further described in varying detail below as follows.

I. Phase Correction

One aspect of the present disclosure provides an improved automaticphase correction approach for MRS signal post-processing. This providescertain perceived benefits versus other approaches previously describedand used, in particular at least one such prior disclosure based onfitting a polynomial, typically a first order or linear, to the phasesequence and then determine the time-zero intercept point which would betaken as the zero-order phase error. The phase correction approach,according to both prior disclosures and the present embodiment, areconfigured to generally operate on the unsuppressed water FID or on theaveraged set of time domain complex free induction decay (FID) frames.

The present phase correction approach is configured to analyze theabsorption spectrum and use an optimization approach to maximize orotherwise increase the real part of the absorption spectrum. Suchoptimization is satisfied by relative enhancements or increases, anddoes not require reaching a true optimum value. This automaticallyperforms a similar function which is otherwise typically required to bedone manually on MRI system consoles or other processors inpost-processing MRS acquired data. The optimization approach accordingto the specific detailed embodiment presently described is based onusing (a) a golden section search, generally described as a techniquefor finding the extremum (minimum or maximum) of a unimodal function bysuccessively narrowing the range of values inside which the extremum isknown to exist (see, e.g., Wikipedia description for some background),and (b) parabolic interpolation, such as available and may be used viathe “MATLAB™” function called “fminbnd,” as may be employed by one ofordinary skill according to the unique applications describedcontextually hereunder. This computes an FFT on the input data withinits iteration loop, though it converges very rapidly. An objectivefunction is to find the phase shift which reduces or minimizes thenumber of FFT bins in the absorption spectrum which have a negativevalue.

Other more complex objective functions may be employed, and have beenevaluated, but generally this approach has been observed to be adequatedespite its simplicity and is not significantly improved upon by morecomplex methods in test cases observed.

For the purpose of further illustration, the algorithm andimplementation are described in the header of one exemplary source codeembodiment provided below.

function [phi, zphi, fftzph, fval, status] =optPhaseAdj(z,phiLo,PhiHi,maxIter ) % Determine the phase angle whichmaximizes the real part of Fourier transform % of the complex vector, z.Search for this value between phiLo and phiHi. The % result will be thephase error used to perform zero-order phase correction of % an MRS FID,z. % Algorithm: Use the MATLAB minimization function fminbnd to find thephase % shift which minimizes the number of fft bins which have anegative value. % The objective function simply phase shifts everysample of the input array by % the trial value of phi, takes the FFT,and returns the count of bins below % zero. Although simplistic, this isappears to provide an optimal (or at least improved) zero-order % phaseadjustment; this suggests that nothing would be gained by computing %something more complex like the sum of in-phase power, though this mayalso be adequate %  Since the phase adjusted array and its FFT arecomputed and available at % convergence, they are returned. % Arguments:% z Complex vector to be phase adjusted in the MRS absorption spectrumsense % phiLo, phiHi Search bounds, in radians for the optimum phaseadjustment % maxIter Maximum number of iterations to allow % Returned: %phi Phase error, in radians, for optimal adjustment % zphi Phaseadjusted complex vector: z with phase adjustment phi applied %fftzph Complex spectrum of phase adjusted vector % fval Value ofobjective function % status Convergence status % 1 Converged tospecified tolerance % 0 Reached maximum iteration count % −1 Terminatedby output function (see MATLAB help) % −2 Inconsistent bounds (phiLo >phiHi) % 20111101 PK % Notes: % 1: Tolerance set to 0.01 - this wasfound to work well on poor SNR data, %  converging in 10 to 15iterations. % 2: maxIter may vary, though start with 25. % 3: MATLABhelp warns that convergence may be slow if solution is close to % search boundary. In any case, if it can be given a narrower bound than−pi %  to +pi, that will help. status = 0; [phi,fval,status] = ...fminbnd(@phaseAdjMin,−pi,pi,optimset(‘TolX’,0.01,‘MaxIter’,maxIter));function fval = phaseAdjMin(ph )  [Th,Rh] = cart2pol(real(z),imag(z)); Th = Th − ph;  [X,Y] = pol2cart(Th,Rh);  zphi = complex(X,Y);  fftzph =fft(zphi);  fval = length(find(real(fftzph) < 0)); end fftzph =fftshift(fftzph); end

II. Frame Editing

Generally speaking, the frame editing process measures characteristicsof each frame (a single time domain Free Induction Decay or “FID”signal). These characteristics are compared to thresholds. Frames thatpass these tests are used to frequency correct previously phasecorrected frames; frames that don't pass are discarded. If aninsufficient number of frames pass the criteria, per an additionalthreshold value that may be chosen and defined, then frame editing onthe whole may be discarded (as may be frequency correction alsodiscarded, if the quality of the retained frames are not consideredenabling for robust frequency correction). These frames are eventuallyused in aggregate resulting spectra to provide signals of sufficientsignal:noise ratio (SNR) to enable quantification of the chemicals ormetabolites represented in the signals. Frame editing removes “lowquality” frames for the purpose of improving the SNR and, ultimately,improving the quantification of the metabolites. Frame editing isapplied to all phase corrected frames prior to frequency correction aspart of the “Signal Quality Analysis” block in FIG. 2. If there aresufficient frames left after editing, then frequency correction isperformed.

According to the present embodiment, and as shown in further flowdiagram of FIG. 3, there are three frame editing modules or criteria:(1) Full Width Half Max (FWHM); (2) Confidence; (3) Frequency Error.Frames with FWHM values greater than a threshold value, confidence lowerthan a threshold value, or Frequency Errors greater than a thresholdvalue, are flagged to be edited or removed from further analysis. In onefurther more detailed embodiment, these three thresholds values are setstatically based upon empirical analysis of a test data set and as aresult of non-real-time (e.g. non-dynamic relative to a given inputsignal) analysis of the dependency of these thresholds on the desiredoutput; high SNR and accurate metabolite measurements. Additionally, asingle set of thresholds is used for all MRS data, independent ofsystem. The thresholds are generally set for intended optimization forthe best expected performance over the expected range of MRS data, basedupon historical data reviewed. However, as may also be elsewhere statedherein, descriptions related to intended optimizations, maximizations,or minimizations of certain parameters of metrics are not absolute andare not required to be fully optimum, maximum, or minimum within theintended scope of the respective disclosure where such language is used.It is contemplated that actual performance may be improved or enhanced,e.g. increased or decreased as a given desire may be, versus otheroptions. Moreover, such performance measures even as to improvements orenhancements may not apply for every case in which the disclosedinvention is put to use. The disclosed embodiments are intended to beapplicable in wide use, and under variable conditions and uniqueconsiderations between cases. Accordingly, the various performanceimprovements or enhancements intended by various aspects of thepresently disclosed embodiments may actually manifest as intended insome cases, but not others.

Full Width Half Max (FWHM) Frame Editing Criteria

The FWHM frame editing embodiment and criteria measurement quantifies aspectral peak. More specifically, an example of FWHM is illustrated byreference to FIG. 4 (which is generally well understood in the art,though shown here as similarly illustrated in the commercial MR scannermanual available from Siemens Corporation, “SIEMENS MR SpectroscopyOperator Manual, Version syngo MR 2002B”). FIG. 4 shows spectralanalysis measurements incorporated into this present embodiment, whichmeasures the FWHM of the water signal, wherein by reference to theFigure:

v=frequency

v₀=frequency of interest

h=peak amplitude

b=FWHM; the peak width at half height

According to one present embodiment, if the water FWHM is greater than athreshold value, then the frame is flagged to be removed (edited) fromthe channel. This editing step is designed to edit out relatively widewater signals, as such width has been observed to vary over the courseof a multi-frame acquisition series. Such wider water signals may atcertain times provide tails which can shroud adjacentchemical/metabolite peak regions of interest, and have been observed torespond unfavorably to water suppression techniques, in particular atcertain times resulting in side lobe artifacts. This frame editingapproach is thus configured in order to prevent such energy in the widewater signals from compromising other metabolite bands. The approach ofthis embodiment looks in a band around the water signal for the peak;finds the two half height points on either side of the peak, andcomputes the width at these points.

One example of this measurement taken across one illustrativemulti-frame acquisition series is shown in FIG. 5. Here an editingcriteria threshold value of about 30 is provided for illustration. Thisexample (and others shown below in further illustration of the presentframe editing embodiments), is illustrative of a single voxelspectroscopy (SVS) acquisition from an intervertebral disc in the lumbarspine of a human subject, acquired via a 3 T MR system and method. Theexemplary threshold value of about 30 has been observed to appropriatelyrepresent a transition between generally acceptable and typicallycompromised spectra in a clinical test population evaluated via similarSVS acquisitions in discs.

Additional comparative examples are shown in FIGS. 6A-D and FIGS. 7A-D,respectively, for two alternative embodiments featured in relation toanother additional test case, also for a 3 T SVS acquisition in a lumbardisc.

More specifically, FIGS. 6A-D show various aspects of one embodimentapplying a FWHM threshold criterion of 50. As shown by the FWHM plotalong the multi-frame SVS acquisition series in FIG. 6A, the FWHMmeasurement becomes regularly compromised after about 80 FID frames intothe acquisition. However, due to the threshold criteria set to 50, mostof these compromised frames are retained—as is shown in the timeintensity plots for the acquisition pre- (FIG. 6B) and post- (FIG. 6C)frequency correction following frame editing according to thisembodiment and FWHM threshold criteria. The resulting post-processedspectrum (including other post-processing modules applied consistentwith other aspects of this disclosure) is shown in FIG. 6D.

For illustrative comparison, FIGS. 7A-D show similar views of the sameaspects for the same SVS acquisition through the post-processing aspectsrevealed, but in this different embodiment a FWHM threshold criteria of20 is applied. The result captures for selective editing/removal manymore of the compromised frames with widened FWHM which exceed thethreshold criteria of 20, as occur beyond the 80^(th) frame, than in theprevious embodiment. This different result is further shown in comparingthe time intensity plots of the frame data taken pre- (FIG. 7B) andpost- (FIG. 7C) frame editing and frequency correction. The fullypost-processed spectrum resulting according to this embodiment, andsimilar additional post-processing as according to the immediatepreceding embodiment, is shown in FIG. 7D. This shows higher peak value,and narrower more coherent line width, for the n-acetyl (NAA) peakregion around 2 ppm along the chemical shift (CS) spectrum, versus theprior embodiment applying the higher FWHM threshold criteria to the sameinput data and otherwise common post-processing approaches betweenembodiments.

Confidence Interval Frame Editing Criteria

The confidence interval is a tool used to signify the reliability of anestimate in a statistical analysis. The interval is the range in asample distribution between which it is expected that the populationvalue will lie, given the particular degree of confidence. Theconfidence interval value (CIV) for a frame reflects the confidence thatthe water spectral line for that frame can be identified fromsurrounding spectral lines. This is done by characterizing the energy(area under the curve) of the water line with respect to the regionwhere the water line is expected.

The CIV measurement starts with determining the FWHM of the spectralpeak (assumed to be water) in the center of a region. The bins (unit ofintegration) within the FWHM of the water spectral line define the waterline energy. The number of bins in the region surrounding and includingthe water signal greater than bins within the FWHM is counted. The CIVis the number of bins within the water line divided by this difference,and ranges between 0 and 1. More specifically, the CIV is calculatedaccording to the following formula:CIV=Water line energy/Total energy in region−Water line energy

As applied according to the present embodiment, frames with relativelylittle energy outside the FWHM of the water line have a higherconfidence than frames with relatively more energy outside this range.

A comparative example demonstrating the beneficial impact of applyingsuch a CIV threshold criteria for frame editing is provided as follows.

FIGS. 8A-D illustrate yet another frame editing embodiment as applied tothe same multi-frame SVS acquisition series from the same intervertebraldisc featured in FIGS. 6A-7D. However, these current figuresillustrating this present embodiment instead reveals that, during themulti-frame acquisition series, the confidence estimate for the watersignal also became compromised. According to this specific embodiment, aCIV threshold of 0.5 was applied for confidence interval based frameediting. More specifically, FIG. 8A shows the confidence of the framesin the series, with significant reduction in the confidence revealed atabout 80 frames into the acquisition (similarly where the FWHMmeasurements became compromised). However, as shown here with the CIVvalue threshold set to 0.5, only a few of these later compromised frameswere also dropped below this threshold value and were edited forremoval. Most of the compromised frames were retained. FIGS. 8B-C showtime-intensity plots for the acquisition series prior to (FIG. 8B) andpost (FIG. 8C) frequency correction following frame editing according tothe confidence interval-based criteria at the 0.5 threshold valuedescribed. The final post-processed spectrum (also incorporating otherpost-processing approaches also herein described following frequencycorrection) is shown in FIG. 8D.

For comparison, the same acquisition series illustrated in FIGS. 8A-D isalso illustrated in similar views in FIGS. 9A-D. However, these figuresillustrate an alternative processing embodiment wherein the confidenceinterval (CIV) threshold is set to 1 (vs. 0.5 according to theimmediately preceding embodiment). As shown in the confidence plot inFIG. 9A, the vast majority of the later frames with compromisedconfidence drop below the 1 CIV threshold after the 80^(th) frame. Theseare thus flagged and edited, with far fewer frames retained, as revealedin the time intensity plot post-frame editing shown in FIG. 9C (vs.prior to frame editing, as shown again in FIG. 9B). The finalpost-processed spectrum according to this embodiment is shown in FIG.9D, and reveals are more coherent NAA peak with narrowed line width inthe CS region of 2, versus the spectrum shown in FIG. 8D for the sameacquisition but processed according to that preceding embodiment withthe different and more forgiving CIV threshold criteria.

Frequency Error Frame Editing Criteria

The frequency error is the absolute value of the difference betweenwhere the water signal peak is located along the chemical shift spectrumand where it should be. If this difference is large then energy from thewater signal can leak into and distort the metabolite peak regions ofinterest, and ultimately compromise robustness of the metabolic spectralrange and corresponding diagnostic interpretation. Leakage is due to thefinite water suppression bandwidth. Frames that exceed +/−maximumfrequency error threshold criteria are flagged for removal from theseries for further processing. Frames that do not exceed the +/−maximumfrequency error threshold are retained and subsequently frequencycorrected to improve coherence (Frequency Correction). As with othercriteria, however, if insufficient frames are retained via thiscriteria, then according to one embodiment all frames may instead beretained and frequency correction is aborted in the overall processingregimen. However, due to the poor signal quality according to thisapplied criterion, the result may be flagged as potentially compromised.

Real-Time Frame Editing Optimization

Another aspect of the present disclosure provides a novel approach foroptimizing frame editing criteria. This embodiment proposes enhancementto the otherwise “static” empirically derived approach described above,by providing an approach that sets the frame editing thresholdsdescribed “dynamically” in real-time based upon optimizing results for agiven input signal. This provides the ability to optimize the frameediting to achieve the highest SNR and best quantification of themetabolites for each input series. Moreover, the optimization criteriamay be, but do not have to be, limited to SNR.

The table below lists the range of thresholds used for optimization. Thenotation that defines the values is of the format; Start:Increment:Stop.The Confidence Threshold, for example takes on 3 values: 0.5, 0.6, and0.7.

The result of these tests would be a 3×5×4×7 matrix space. Theoptimization search would then be to find the set of thresholds thatmaximized the desired metric, SNR for example.

III. Frequency Correction

Another aspect of the present disclosure is an improved approach forfrequency error correction in MRS spectral signal post-processing. Thefollowing describes the processing flow and functions performed by theMRS post-processor to effect robust frequency error correction accordingto this present embodiment.

A primary purpose of frequency error correction is to reverse frequencyshifts which occur on a frame-to-frame basis, such as for example as mayoccur due to subject movement or to Bo variations resulting from subjectrespiration. Static frequency errors also occur if the MRI system wasnot set to correctly center the water signal to zero at the start of theseries acquisition. A primary objective of post-processing frequencycorrection is that the residual water signal be positioned at “d.c.” (or0 Hz) in the spectrum for each frame. When the appropriately centeredframes are combined by averaging, the metabolite signal components willbe aligned in frequency and phase and will therefore sum coherently—or“coherent averaging.” This coherent summation, post frequencycorrection, will often result in increased metabolite peak height andnarrower linewidths. This approach minimizes signal loss due tonon-coherent averaging and linewidth spread due to frequency errorartifact. As the noise component is more random (relative to the signalcomponent in the frame spectra), the frequency alignment increases thesignal without similar increase in noise, thus also having the result ofincreased signal:noise ratio (SNR).

The approach of the present embodiment is described in further detail asfollows. First, the MRS multi-frame acquisition series is prepared forfrequency correction, such as for example employing one or more (or all)of the following: phase correction, channel selection, frame selection,zero-padding, and apodization. Frequency error is then estimated, todetermine the frequency error for each frame. Frequency correction isthen conducted to apply a correction to negate the frequency error.Coherent averaging is then conducted of the complex frequency-correctedframes, thus computing the final composite spectrum as coherentlyadjusted according to this novel method.

Frequency Shift Estimation

More detailed frequency estimation embodiments are described as follows.One embodiment for estimating the frequency error operates on thediscrete Fourier transform absorption spectrum of each frame todetermine the frequency bin representing the peak of the residual watersignal. The analysis is limited to a parameter-specified range aroundthe zero frequency point or “center bin.” For one exemplary embodimentconfigured to work with MRS data acquired via a 3 T MR system, thisrange corresponds to plus and minus about 50 Hz.

In addition to determining the peak frequency of the residual watersignal, several signal quality and reliability estimations are madewhich are used in the frame selection process to determine if thefrequency error estimate is of sufficient reliability to be applied.These measures include a numerical confidence and an estimate of thelinewidth defined as the full width at half amplitude, in Hz—both of thewater signal observed. The result of the frequency error estimationprocess is a vector of error estimates, in Hz, corresponding to each FIDor MRS frame.

The confidence measure according to the foregoing is computed accordingto another more detailed illustrative embodiment as follows. Thediscrete amplitude spectrum is analyzed in the range of the center-tunedfrequency plus and minus a threshold value, such as for example about 40Hz (which may be appropriate in the case of data acquired via a 3 Tsystem, and about half that for example with respect to data acquiredvia a 1.5 T system). The highest peak is located within the range, andits width at the half-amplitude point is determined. Next, the totalspectral width of all parts of the spectrum which exceed thehalf-amplitude point of the highest peak is determined. The confidenceestimate is made by taking the ratio of the spectral width of thegreatest peak divided by the total spectral width which exceeds thethreshold. If there is only a single peak above the threshold, theconfidence estimate will be 1.0. As a number of peaks or spectralcomponents which could be confused with the greatest one increases, thenthe estimate will be reduced accordingly toward approaching 0.0. Thisprovides a simple and robust estimate of the randomness or dispersal ofenergy in the vicinity of the water peak. Like an entropy measure, ithas the desirable characteristic that its performance is generallyinvariant with amplitude.

An alternative approach believed to offer potentially improvedperformance, in particular in cases of relatively lower water signals inspectra, is based on spectral correlation. According to this approach, acomplex cross-correlation sequence is computed of the complex apodizeddiscrete fourier transform (DFT) of each FID and a reference spectrum.The reference spectrum can be that of one of the individual FIDs, or theaveraged spectrum of multiple FIDs. The process uses the informationinherent in all the spectral components which are common and coherentbetween the test spectrum and the reference to determine the spectralalignment of the two spectra. This is indicated by the location of thepeak of a cross-correlation coefficient. One advantage contemplated forthis embodiment, such as in some cases versus the spectral analysisapproach described above, is that in the case of closely locatedcoherent spectral artifact signals, which has been observed in someacquired spectral examples, this approach uses this artifact signal toadvantage rather than it being disruptive as with spectral analysis.More specifically, the spectral analysis approach above may confuse theartifact with the water signal, decreasing confidence level. Thisalternative approach benefits from all correlating signals. Thecorrelation method may also offer greater sensitivity than the spectralanalysis method in certain cases, as the preceding approach exploitsonly the residual water signal. For example, given sufficiently strongmetabolite signals, the current embodiment of spectral correlation couldbe used with complete water suppression. It is to be appreciated,however, that in relatively low metabolite SNR examples (such as smallvoxels), such metabolite signals may not be sufficient to provide thisadvantage in some circumstances—such that the water-based spectralanalysis approach may be equivalent (or even improved in some cases). Anadditional advantage of the present embodiment, however, is that thecorrelation peaks are better behaved than spectral peaks of residualwater in the presence of artifact, making peak analysis more precise.The result of this process is provided in the same format as for thespectral estimate approach.

Frequency Shift Correction

Given the frequency error estimates for each frame, such as determinedby the frequency error estimate approaches described above, thecorrection is then applied to the data. This is done according to onepresent embodiment in the time domain.

According to another more detailed embodiment, this is done bymultiplying complex FID data by a linear complex phase function, ofunity amplitude, representing the negative of the frequency error. Thisis a function for which the phase increments by (2·π·−ferrHz/Fs) foreach sample, where ferrrHz is the frequency error in Hz, and Fs is thesample rate in samples per second. This has the effect of rotating thespectrum such that the measured water peak will be shifted to the zeroor d.c. location. This linear function also has effect of performing afirst-order phase correction across the spectrum. This frequencycorrection operation is applied to the frames of each acquisition (coil)channel which are to be combined.

As previously described, it is further contemplated that frequency errorcorrection can also be conducted in the frequency domain—simply byshifting the spectra by the appropriate bins to coherently align them(e.g. via the water peak, per estimate of the bins represented by theshift error from the center point for each FID).

IV. Frame Averaging

Frame averaging is performed on coherently aligned FID frames, per theforegoing approaches. This is done according to another presentembodiment in the time domain, and may be accomplished by simplyaveraging the selected complex frame data for each FID. This isperformed for the FIDs collected from each of the channels which are tobe combined. It is further contemplated according to still furtherembodiments that this could also be well performed in the frequencydomain, though the time domain approach has been observed to be moreefficient in many cases.

V. Phase Grouping

For example, another such aspect comprises processing MRS data from anacquisition series conducted across first and second phase groupscomprising first and second sets of acquisition frames acquired at firstand second different respective phases along a phase cycle. Further tothis aspect, the system and method processes at least in part the firstand second phase cycle groups separately. According to one mode of thisaspect, the separate phase group processing comprises at least one ofthe other processing aspects, modes, embodiments, variations, andfeatures elsewhere herein described, such as for example one or more offrequency error estimation, frame editing, frequency correction, andphase correction. According to another mode, the first and second set ofacquisition frames for each phase group are combined into one phasegroup combined result. According to another mode, the first and secondphase groups are combined to provide an averaged spectral resultfollowing the separate processing. In another mode, the first and secondsets of each phase group are combined (e.g. averaged), and then thecombined phase groups are combined (e.g. averaged). Another modecomprises reducing or removing artifact separately between the phasegroups. In another mode, the artifact is removed by combining the phasegroups after separate processing and in-group frame combining. Inanother mode, frame editing is performed on each phase group, and adifferent number of frames are retained (or conversely filtered out) ineach group prior to combining in-group frames and the groups together.

Another aspect provides an MRS acquisition data set in separate phasecycle groups, and a process which is performed on each phase cycle groupindependently of the other. According to one mode of this aspect, thephase of certain excitation pulses in the pulse sequence are changedsuch that unwanted coherences will have different phases when the signalis acquired and while desired coherences will have the same phase, Thisis so that when the frames are combined, the artifact signals cancel andintended signals combine constructively. It is to be appreciated thatsuch phase step cycling could be transparent to post-processors andstill achieve sufficiently adequate results in many circumstances. It isalso appreciated, however, that certain benefits may also be provided byintegrating phase cycle grouping into post processing, such as forexample to mitigate degradation of results due to artifact. According toone example embodiment for illustration, a processing approach usingframe editing may result in fewer frames in one phase cycle group thanin the other, versus matching them as typically intended for equallyweighted averaging. By processing the frames as one overall groupwithout regard to phase cycle grouping, the group with more frames wouldhave a higher weighted contribution to the averaged result. This mayresult in a phase bias in the result with regards to unwanted signalqualities (e.g. artifact) intended to be removed by the very purpose ofthe phase cycling. By treating the phase cycle groups separately, andaveraging them separately, their respectively averaged resultscontribute equal weighting to then thereafter averaging those interimresults together. In addition, other processing operations such asspectral correlation may not perform as well when attempting tocross-correlate FID frames of separate phase cycle groups (due to theirrespective phase differences), vs. within the same phase cycle groupsfirst.

VI. Channel Combining

The MRS acquisitions are typically collected on multiple channelsrepresenting the signals acquired from a multi-coil antenna detectorarray. Zero-order phase alignment and frequency correction (with orwithout frame editing) will have been performed for each channel, thoughchannel selection may have been employed to reduce the channels from allwhich acquired signals to only those determined to be the strongest (perother embodiments contemplated). The remaining step is to combine thecomposite FIDs from each channel retained. This may be performedaccording to one approach by merely combining them into one compositeaverage. This may suitably approach optimal results when combiningchannels of similar SNR. Alternatively, this may be done using amaximal-ratio combining approach, such as often employed formulti-channel radio frequency signals in communications technology.According to one illustrative detailed embodiment, consideredparticularly beneficial, the channel weighting factors used are based onthe total power in each retained channel after apodization.

VII. Mirrored Downfield Artifact Correction in Upfield Metabolite Range

Another aspect of the present disclosure relates to signal componentsinduced MRS spectra downfield of the water signal in the spectra. Thisis described by way of various particular embodiments below, includingby reference to certain Examples 1-4, and corresponding Figures byreference thereto. Such induced downfield signals, when present, havebeen observed to co-exist with similarly located signals mirrored in theupfield side of the water signal. When present, the locations of certainof these downfield signals have been observed with appreciableregularity to roughly mirror the locations of metabolite ranges ofinterest upfield of water. This is particularly observed at downfieldregions roughly mirroring peak locations opposite the water reference orcenter line corresponding with at least one of n-acetyl (NAA) (e.g.1.9-2.3 mirrored location), lactic acid/alanine (LAAL) (e.g. 1-1.5), andlipid (e.g. 0.5-1.5) chemical regions along the spectrum. When thesedownfield signals are generated, they also correspond quite regularlywith similarly shaped peaks generated at these respectively mirroredlocations.

An exception is found in relatively discrete down field peaks that areoften generated more closely downfield of water, e.g. 4.7-6.1 ppm in theCS spectrum. These induced peaks closely downfield of water have notbeen observed to correspond with mirrored peak signals on thecorresponding closely upfield side of water (e.g. 3.3-4.7 ppm). It isbelieved that this closely downfield peak, when present, may correspondwith one or more MRS chemical signatures unique to that location (e.g.uptake of certain pain medications may manifest in such peaks).

It is believed, based on scientific theory and also certain supportingobservations (including via observations in clinical test data), thatthese downfield signals and corresponding mirrored upfield signalsaround water (e.g. beyond about 6 or 6.1 ppm, and below about 3.3 ppm,respectively, in CS spectrum) may comprise side lobes induced byapplying water suppression during an applied MRS pulse sequence to, andrelated MRS signal acquisition from, a prescribed tissue region ofinterest in a patient. This would support the roughly symmetricallymirrored and similarly shaped signals on each side of water/center line.Because no appreciable chemical or metabolite signals are expected to begenerated in this downfield range, the signals observed there areconsidered to be the downfield components of such artifact. Often thecorresponding upfield peaks generated in the metabolite regions ofinterest have more power than these mirrored downfield peaks. However,these upfield regions are expected to comprise real, discrete chemicalpeaks known to correspond with the chemical constituents noted (when intissue)—and so are expected to be additive between real signal (whenpresent) and the reflected artifact noted.

According to the foregoing, this present aspect of the disclosure usesspectral peak artifacts induced in the downfield CS spectrum as a basisby which to adjust the upfield spectrum in a corresponding range roughlymirrored around the water peak location (about 4.7). This is done toremove suspected contributions of similarly mirrored artifact componentsfrom real induced chemical signal peak components that may be collocatedin this upfield range. This is generally done by subtracting the poweralong the downfield range from the corresponding mirrored upfield range,as further described below. According to one mode, the downfield rangecontemplated is from about 6 or 6.1 ppm and beyond, and the upfieldrange is from about 3.3 ppm and below, for this adjustment.

Accordingly, one current embodiment captures the downfield signal,inverts it around an axis (e.g. water) to align it against the upfieldspectrum, and subtracts the inverted downfield spectrum from the upfieldspectrum to provide a resulting “downfield adjusted” spectrum in theupfield range. It is further recognized that the “mirroring” effect ofthe downfield artifact signals may not be exactly perfectly symmetricalin all cases, and so simple inversion around the water axis may notresult in exact corresponding alignment between the inverted downfieldartifact peaks and the upfield peaks. Among other reasons, side lobesarising from water suppression anomalies would not necessarily beperfectly symmetrical, especially via water suppression conducted via afixed range around 4.7 ppm in the setting of an actual water signal thatmay not exist exactly at 4.7 (or vary about it between frames due tofrequency shift and/or other anomalous influences). Accordingly, thealignment of the reflected downfield spectrum with the upfield spectrummay be further adjusted to better align the peaks—to the extentresulting in enhanced correlation, and more robust expected spectralresults arising from the artificial subtraction step imparted on theoriginally induced spectrum.

The foregoing general description is further developed in finer detailaccording to the additional exemplary embodiments as follows:

More specifically, downfield (generally referring to chemical shiftvalues greater than water; nominally 4.7 PPM) artifact spectral peaksignals believed to be related to water suppression has been observed insome of the clinically acquired MRS data and observed to also correlateto similar artifact signals quasi-symmetrically located on opposite “upfield” side of water in range of interest for metabolite signal regions.In particular, downfield peaks have been observed quasi-symmetricallyaround water to correspond with proteoglycan (PG) and lacticacid/alanine (LAAL) related regions of diagnostic interest (among otherregions, as elsewhere noted herein). The purpose of the followingexamples is to uniquely reflect observed downfield artifact intoup-field metabolite regions, and artificially apply innovative solutionsto reverse and correct for this in spectral presentation and/or analysisfor diagnostic purposes.

Adjustment of the up-field artifacts by the down-field spectralcomponents, believed to be symmetrically reflected as embedded in theupfield spectral range, is conducted according to the following steps,with the adjustment process illustrated using plots for four exemplarydisc studies further described below:

1. Coarse Alignment of the spectrum to DC

2. Fine Alignment of the up and down field spectral components in aspecific chemical shift range

3. Baseline Correction (an optional step, and may be done prior to orafter adjustment).

4. Adjustment: subtraction of the aligned down field components from theup field components in a specific chemical shift range of interest.

i. Coarse Alignment

Coarse alignment of upfield and downfield spectra is performed accordingto one embodiment as follows. The location of the water peak isdetermined. This is typically the maximum peak value in the spectrum andshould be located in a DC-centered spectrum at or at least very near bin1025 (e.g. when using 2048 point FFTs). All plots shown are in terms ofChemical Shift in PPM. Water corresponds, roughly, to a chemical shiftof 4.707. If the peak is not at DC then the entire spectrum iscircularly rotated so that the peak is located at DC. FIG. 10 shows achemical shift spectrum and its mirror, using water as the line ofsymmetry. This example shows that the computed water peak is one bin tothe left of water. The mirror is then one bin to the right of water.FIG. 10 shows the input spectrum representing an initial stage ofprocessing a first Example 1 spectrum. The input spectrum shown ispost-processed according to various embodiments herein shown anddescribed, but without baseline correction and pre-downfield artifactcorrection according to the current embodiments.

FIG. 11 is the same spectrum at same stage of processing shown in FIG.10, but showing a broader scale of the spectrum with the water peak anddownfield spectrum to its left in the graph, and with salient side lobes(which may be for example representative of water suppression artifacts)marked to illustrate that the water peak is not always at DC, and alsothat the water suppression artifacts are not always exactly symmetricabout water. It is noted referring to FIGS. 10 and 11 that the signs ofthe values are negative due to the way data is plotted and theleft-to-right decrease of the Chemical Shift axis. Nonetheless, spectralrange to left of water is considered “downfield” whereas spectral rangeto right of water peak is considered “upfield” for purpose of describingthese current embodiments of this disclosure.

FIGS. 12A-B show further progression of Example 1, in sequentiallyincreasing resolution (FIG. 12B showing mainly the metabolite range ofinterest). This illustrates the same input spectrum and its mirror, andillustrates the 1 bin of absolute error of both signals relative to DC.

FIGS. 13A-B show the spectrum circularly rotated one bin so that thecomputed water bin (bin with the maximum value) aligns with DC (ChemicalShift of 4.707 ppm).

ii. Fine Alignment

Next the original and mirrored spectra are more finely aligned. This isbecause the location of down-field and up-field spectral components arenot always exactly symmetric about water, though it is believed thatwhere significant peak correlation results from fine tuning of thealignment, such step is supported as reasonable. Such fine alignment isdone according to the particular present embodiment by correlating thespectral components in the chemical shift range [3.2-0.5]. FIGS. 14A-Bshow the fine aligned, full and zoomed spectra, respectively.Specifically noted by reference to that figure, the edges of the watersuppression artifacts are in alignment. In this illustrative example ofFIGS. 14A-B, the spectral components are aligned when the mirroredspectra is shifted 44 bins to the right. Only the mirrored spectrum isshifted so that we do not change the location of the up-field componentsas they are correctly located and aligned to water.

iii. Baseline Correction

Baseline correction is an optional step, taken in the present embodimentprior to adjustment of the up-field spectral components via reflecteddownfield signal subtraction. The baseline correction algorithm operatesonly on the metabolite range of interest up field of water. A polynomialis used to fit data in the metabolite range. This method will notgenerally work on the entire spectra, though only ranges of interestneed be worked upon in any event. Baseline correction is done separatelyon the up-field and mirrored down-field spectra (mirrored to the samemetabolic range). FIG. 15 shows the baseline corrected spectra (onlydone on up-field spectral components) and FIG. 16 shows the baselinecorrection of the mirrored down-field spectra.

iv. Subtract Downfield from Up Field Spectral Components

At this point the down-field spectral components are aligned with theup-field components and we have baseline corrected versions of both theup and down field spectra. Adjustment is accomplished by subtractingmirrored downfield spectra from the up field except in the up-fieldrange corresponding to the down-field component in the Chemical Shiftrange [from about 4.7 to about 6 or 6.1] PPM. What follow are plots ofthe spectral components and correction using the baseline andnon-baseline corrected spectra.

a. With Baseline Correction

FIG. 17 shows the baseline corrected up-field spectra and the baselinecorrected mirrored down-field spectra after coarse and fine alignment.FIG. 18 shows the adjusted up-field spectra.

b. Without Baseline Correction

FIG. 19 shows the up-field spectra and the mirrored down-field spectraafter coarse and file alignment, but without baseline correcting theserespective spectra prior to adjustment. FIG. 20A-B show the adjustedup-field spectra after applying the downfield subtraction, which isshown pre- (FIG. 20A) and post- (FIG. 20B) baseline correction followingthe adjustment, respectively.

Further examples are similarly illustrated in FIGS. 21A-25C, 26A-30C,and 31A-35C, respectively, for three additional clinical test casesinvolving 3 T acquisitions from other discs at other lumbar levels inother subject exams.

Spectral results are evaluated by plotting the adjusted up fieldspectrum in the range of metabolic interest PPM, and generally limitingthe Y axis values to provide optimal observation between (e.g. +10%) ofthe maximum and minimum in the defined range. Comparison is then madevs. final post-processed spectra without the downfield subtraction tounderstand possible contribution of the artifact being evaluated, andgeneral spectral quality analysis (e.g. SNR). Metrics, such ascomparison of resulting SNR of signal region of interest pre- vs.post-adjustment, or correlation measurement of the reflected aligneddownfield signals with the upfield signals, may also be used accordingto further embodiments to then determine whether the adjusted result, ornon-adjusted result, should be used.

It is also appreciated that measurement data of downfield artifactsignal may be subtracted from upfield metabolic ranges of interest, vs.performing spectral adjustments. This may relate to the same alignmentapproaches as above for spectral adjustment, but simply not change theshape of the spectrum but rather simply adjust the spectral measurementsto account for estimated downfield reflected artifact contribution.Among other benefits of this approach, this does not require a “step”change in the resulting spectrum at the edge(s) where subtraction isused versus not used.

This data subtraction approach may be conducted as follows. Similarrange “alignment” approaches may be used as with spectral subtractionadjustment approach above. However, then this is merely used todetermine ranges of up-field and reflected downfield artifactmeasurements to take for subtraction adjustment. This may involvemeasuring peak and/or areas under the curve (AUC) for metabolic regionsof interest, and reflected similar ranges on downfield side of water(e.g. other side of 4.7, or as reflected and aligned, either coarse orfinely). Then the values for downfield ranges are subtracted fromaligned upfield ranges, providing “corrected” measurements. This mayalso be evaluated versus upfield ranges measured prior to correction todetermine artifact contribution.

The robust basis and approach for reflected downfield artifactcorrection in upfield signal ranges of MRS spectra is supported by wayof this description and illustrative examples shown. It is alsoappreciated that this may be further developed through refinements overlarger sample sizes of data. In particular, custom variations may beemployed to address certain particular types of reflected artifactsignals observed, and establish criteria by which to either usedifferent approaches, and/or whether to use an adjusted spectrum (ormeasurement) or not in particular cases.

v. Validation of Mirroring Function

Subtracting reflected downfield signals from upfield metabolite range ofinterest is done according to the present embodiments principally basedupon the presumed principles that that downfield signal is in factartifact which is mirrored into the upfield range at the place wherereflected downfield is aligned for subtraction. Accordingly, ifmirroring of the input spectrum for downfield correction is notappropriately aligned to where real artifact is, the subtraction willinappropriately adjust the upfield signal by newly induced artifact bythe inappropriately conducted subtraction. Accordingly, the mirroredreflection approach can be validated for proper alignment as intended,which is a useful approach to take as mis-alignment could confoundresults as stated. To validate this function, a further aspect of thepresent disclosure provides validation approach that comprises a simpletest vector that is submitted to the operation, for its mirror image tobe examined. The test vector is designed to simulate an input signal,and which may be for example the real part of a complex spectrum.

One particular example of such a validation is provided herein asfollows. For simplicity of illustration, this example simulates a16-point real spectrum. With this signal, DC is at bin 9 (16/2+1) andthere are 7 (16/2−1) bins upfield from DC and 8 (16/2) bins downfieldfrom DC. As there is not the same number of upfield and downfield bins;there is one “odd” bin, the furthest downfield, without a correspondingupfield bin. One issue thus presented by this approach is how to handlethe “odd” bin. If the input vector were truly mirrored, then themirrored vector would have 7 bins downfield and 8 bins upfield of DC.This would put DC at bin 8. The signals would be mirrored but the DClocations would be different. Alignment of the DC points could be forcedby prefixing the mirrored sequence with one bin with amplitude zero, butthen the test and its mirror would be different lengths. The otheroption, and the method implemented, is to mirror the entire spectrumexcept the “odd” bin. The leaves the DC bins aligned and the vectors arethe same length.

FIG. 27 shows the discrete sampled test vector (“Input”) in dashed lineand its mirror (“tupnI”) in dotted line. The signal is zero except at DC(amplitude=5), DC+5 (amplitude=1), DC−3 (amplitude=2), and the “odd” binDC−8 (amplitude=3). The mirrored test vector shows that the mirroringfunction is correct: The “odd” bin has not moved; the DC bins arealigned, and the upfield and downfield components are mirrored about DC.

FIG. 28 shows the results of a test around bin 12, with center findingfunction validated.

VIII. MRS Pulse Sequence Protocol—Water Suppression & Shimming

Another aspect of the present disclosure provides certain enhancementsto an MRS pulse sequence protocol for optimally acquiring robust spectrain clinical applications, and considered to be particularly beneficialin single voxel spectroscopy (SVS) (though also applicable tomulti-voxel), and in particular to MRS acquisitions in relatively smallvoxels, and still more particularly in musculoskeletal joints. This isespecially beneficial with respect to MRS performed in intervertebraldiscs and furthermore in particular in such discs of discogenic low backpain patients where MRS may be deployed for Dx purposes to help identifypainful versus non-painful discs.

Water varies dramatically in intervertebral discs, in particular betweennormal healthy discs and degenerated discs. Accordingly, in discogeniclow back pain patients in particular, a wide range of water signal powermay be found, including between discs of the same patients. As water isused in certain MRS analysis, and processing, this water variancebetween discs, and water suppression approaches that may be uniformlyappropriate, can be an issue. More specifically, in highly hydrateddiscs, such as normal discs, water content is obviously high. No watersuppression used in MRS pulse sequence studies of such discs may resultin water signal that is too strong and may overshadow the metaboliteranges of interest via “tails” of the water signal. In highly desiccated(dehydrated) degenerated discs, water content is obviously low. Watersuppression may remove any water signal that may exist, thereby to thedetriment of any processing or analysis approaches relying on watersignal.

Accordingly, one aspect of this disclosure provides an MRS pulsesequence exam configuration and protocol that comprises a degree ofwater suppression that is set to “moderate” suppression (e.g. betweenfull suppression and no suppression). This may be achieved, for example,for a third flip angle in a CHESS/PRESS MRS pulse sequence of about 60to about 120 degrees, and in particular about 85+/−10 degrees. It is tobe appreciated though that the appropriate exact settings may varybetween MR system and from different vendors, per differentconsiderations in their overall operation and also specific watersuppressions approaches employed. However, in general, such moderatesuppression has been observed to be appropriate in the most cases.Nonetheless, even this has been observed to be inappropriate in othercases, as either too strong or too weak (per the situations notedabove).

According to another aspect, water suppression is controlled to meet agiven case. In one embodiment, a disc is evaluated as to its degree ofhydration. This may be done for example via review of T2-weightedimaging, T2-mapping, or T1rho imaging and/or mapping, of the disc. Basedupon the degree of hydration, a degree of water suppression is selected.For highly hydrated discs, more aggressive suppression is selected, andfor less hydrated discs lower suppression is selected and set, for theMRS exam configuration and protocol customized to that disc. This may bedone automatically, such as via a look up table correlating hydrationmetric to water suppression level, or manually (either by the sameprocess, or otherwise).

Another issue observed with water suppression is the range upon which itis applied. This may be for example 50 Hz (e.g. in case of 3 T), or 25Hz (e.g. in the case of 1.5 T). While such range is typically providedin default setting around an expected “zero” center line for water(about 4.7 ppm in the CS spectrum), this does not always function asintended. In particular, water side lobes may result outside of thatrange post-suppression—giving rise to issues noted above with respect tosuch lobes providing additive artifact in other metabolite ranges ofinterest outside of the water suppression band. This is in particularthe case where there is a wide water line resulting from an inadequateshim. A combination of a broad shim, and narrow water suppression band,has been observed to produce compromised spectra believed to be due tothis phenomena. Many MRS pulse sequence protocols are configured andconducted using only automatic shimming. However, this is based uponcertain defaults that are not always sufficient, in particular in thehigh susceptibility area of the lumbar spine for example. Observation ofthis in the clinical environment has been further observed to result inthe issues noted above.

Another detrimental effect of the relatively narrow water suppressionpulse is that it can introduce a bias error in the frequency shiftestimator. The spectral shape of the water suppression pulse istypically Gaussian; the spectral shape of the residual water shape issimilar. When the water signal falls on the slope of the watersuppression signal, the spectral shape of the water is distorted and itsapparent peak, as sensed by the frequency shift estimator, will not beits true energy peak. This effect is mitigated by the more gentle slopeof a wider water suppression pulse.

According to another aspect of this disclosure, an MRS pulse sequenceacquisition configuration and protocol provides for a manual shim, whichmay be a default approach for the pulse sequence or override in settingswhere autoshim is insufficient. This has been observed to overcome theshortcomings of the default autoshim in many lumbar disc acquisitions,providing more robust results and often avoiding the interrelated issuesnoted with respect to wider water peaks and a fixed water suppressionband.

Another aspect of this disclosure is an MRS pulse sequence acquisitionconfiguration and protocol which provides for wider water suppression.As opposed to the typical bands noted above, increasing these isbelieved, and has been observed in certain clinical cases, to provideimproved results with respect to reducing water suppression side lobeartifacts, though while not compromising metabolite regions of interest.According to one embodiment, the water suppression band is set to morethan about 50 Hz (for 3 T), and more than about 25 Hz (for 1.5 T). Inone further variation of this embodiment, the band is set to no morethan about 100 Hz (for 3 T) and no more than about 50 Hz (for 1.5 T). Instill a further variation, it is set to about 75 Hz (for 3 T) or about37.5 Hz (for 1.5 T). These values are considered reached within about+/−10 Hz. Specifically, settings of 75 Hz and 100 Hz are both expected(and have both been observed in limited cases tested) to provide robustresulting spectra for highly hydrated discs with strong water signal,and without compromising spectral ranges of interest includingpreserving robust carbohydrate signal (which is particularly exposed ascloser to water, e.g. up to about 3.5 ppm), and without appreciable sidelobe artifact of note.

Notwithstanding these ranges as potential limits, however, even theseranges can expand in certain circumstances. For example, wheremetabolites further outside of water are not necessary for a particularmeasurement purpose, the bandwidth can expand to ensure no waterartifact while still preserving other spectral peaks further away. Forexample, if NAA (e.g. proteoglycan) and/or metabolites to its right(e.g. hypoxia metabolites, e.g. lactic acid and/or alanine, or lipids,etc.) are of interest, and other data between the NAA peak region(around 2 ppm) and water (around 4.7 ppm) are not needed, watersuppression bandwidth can extend much wider than these limits up to aspan that would compromise the NAA signal region.

While this aspect of increased water suppression band is believed toprovide a robust solution, in particular for disc MRS (and especiallylumbar), a further embodiment of the present disclosure provides for anadjustable water suppression band in custom manner to accommodate agiven test case. According to one further more specific embodiment, theband is adjusted based upon the value of at least one parameter of thewater signal, such as for example peak amplitude and/or FWHM of the shimresult. According to another more detailed embodiment, a default watersuppression band setting is used for an initial acquisition sufficientto evaluate the result and determine an appropriate adjustment for thesetting. This may relate to a “trial and error” approach around certaindesired result. In still another more particular embodiment, thedownstream field of the spectrum is analyzed for peaks considered to beside lobe artifacts of suppression. The suppression is adjusted to avalue to reduce this artifact. Other metrics, such as upfield effects,may also be used (either alone or in combination).

IX. MRS Diagnostic Processing—Intervertebral Discs

Other aspects of the present disclosure provide certain new and improvedsystems and methods for diagnostic interpretive processing of MRSspectral information from MRS spectra that are post-processed (such asfor example according to one or more, or all, of the various signalprocessing aspects disclosed herein) from MRS pulse sequenceapplications and induced signal acquisitions in and from intervertebraldiscs, respectively.

One such aspect measures one or more parameters (e.g. peak height, areaunder the curve or AUC, FWHM, etc.) for n-acetyl (NAA) and carbohydrate(CA) peak regions, and computes a ratio therebetween. According to onemode of this aspect, this ratio is combined with at least one additionalratio comprising such measurements for at least one of these regions anda second measurement comprising the lactic acid region (LA). In oneembodiment of this mode, the second measurement comprises the combinedlactic acid and alanine (AL) (together “LAAL”) region. The multiplerelative ratios between these several peak regions are provided inrelative context, which is expected to provide diagnostic value incertain cases, such as for example in diagnosing degenerative discdisease and/or discogenic low back pain based there upon and via theuniquely induced MRS pulse sequence signals acquired, post-processed,and measured (such as according to the various other embodiments hereindisclosed).

Another aspect of this disclosure calculates a “peakiness” value for thecarbohydrate region, which is believed to relate to degree of complexcarbohydrate breakdown into constituent component molecules. This isdone according to one mode by curve fitting the region and calculating acorrelation coefficient therefrom. This calculated measurement may alsobe used, either alone or in combination with other metrics, in a manneruseful for assisting in diagnosing disc disease or pain.

It is also appreciated that each such measured spectral feature mayprovide valuable information on its own and without necessarily beingrequired to be taken in combination with other features (such as per theillustrative embodiments described herein). Moreover, other chemicalregions may be of diagnostic interest in certain applications, alsoeither alone or in various combinations, as would be apparent to one ofordinary skill.

IX. Phase Group Processing & Combinations of Embodiments

The following reference is herein incorporated in its entirety byreference thereto: Bolan, Patrick J., et al., “Measurement andCorrection of Respiration-Induced Bo Variations in Breast ¹H MRS at 4Tesla,” Magnetic Resonance in Medicine 52:000-000 (2004).

Additional embodiments including phase step cycling of MRS serial framedata acquisition and related processing, and certain combinations ofvarious embodiments herein described, are further disclosed as follows.An MRS system is configured to perform an MRS spectral acquisition withphase changes stepped along a cycle corresponding to a sub-set ofsequential FID frames. This is typically repeated over the course of theacquisition, such that the total number F of FID frames in theacquisition divided by the number of steps S in the cycle will typicallybe equal to a whole number. A post-processor is then configured topost-process FID frames comprising similar phase steps within “groups”,though the result of which processing may be later combined for furtherprocessing.

MRS systems will frequently employ phase cycling, wherein serialacquisition frames are stepped through different phases. These aretypically done by even numbers of “phase groups” across the acquisitionseries. Certain additional present embodiments that are consideredadvantageous, although not necessary in all cases, integrate frequencyshift estimation (& subsequent correction) and frame selection withphase group processing. The broad scope of the present disclosurecontemplates a variety of approaches for this beneficial processingmethod (and related systems and processors), as would be apparent to oneof ordinary skill.

One particular example embodiment however integrates a spectralcross-correlation based frequency estimator with phase group processingand frame selection. These are integrated according to this embodimentfor two primary reasons.

(1) The frames from each phase group contain artifact with differentphasing. This artifact can introduce a phase-group dependent bias in ourfrequency shift estimator. For example, if one is +2 Hz and the other is−2 Hz, and the correction is applied with this bias, then when theframes are combined, the composite spectrum will be the combination ofthe two shifted spectra and split peaks mat result. This is avoided byfrequency adjusting the groups separately, then performing an absolutefrequency correction on each before combining.

(2) Frame selection can adversely interact with phase grouping. Thebasis of achieving artifact cancellation via phase cycling is that anequal number of frames from each phase group be combined. Frameselection can upset that balance by upsetting that equality andintroducing a bias. This is avoided by having the last step be formingthe mean of each phase group, and then the mean of those means. A testmust be performed to verify that each group is adequately represented,otherwise, phase groups are not processed individually and all theframes are treated as if from a single phase group.

The following describes integrating a spectral cross-correlation basedfrequency estimator with phase group processing and frame selection.These are integrated for two primary reasons.

The frames from each phase group contain artifact with differentphasing. This artifact can introduce a phase-group dependent bias in ourfrequency shift estimator. For example, if one is +2 Hz and the other is−2 Hz, and the correction is applied with this bias, then when theframes are combined, the composite spectrum will be the combination ofthe two shifted spectra and split peaks mat result. This is avoided byfrequency adjusting the groups separately, then performing an absolutefrequency correction on each before combining.

FIG. 38 shows a spectrum from a lumbar disc acquisition in overlaybetween multiple different aspects of the acquired data as follows. Thisis an example of the differences in the chemical shift spectra from aseries of acquisitions with 2-step phase cycling. The dotted line is theaverage of 128 frames from one phase group, and the dashed line is theaverage of 128 frames from the other phase group. The portions of thespectra where these traces are of opposite polarity represent theartifact components that phase cycling attempts to mitigate. The solidline shows the average of these two sep″/showing the cancellation,albeit imperfect, of these artifact signals.

FIG. 39 shows a flow diagram of an automated method to perform frequencyshift estimation and correction, and frame selection in the context ofphase group processing. It uses the spectral cross-correlation method offrequency shift estimation but could be based on the spectral peakmethod just as well. The first step is to form reference spectra foreach phase group. For the second step, all the frames from thecorresponding phase group will be aligned in frequency by the spectralcross-correlation method. The third associates frame selection and phasegroup membership to form the mean of the frequency-corrected andframe-selected frames for each phase group, and then the mean of thesemeans. This is performed in the time domain (TD) and then transformed,with appropriate apodization to the frequency domain (FD) to form thefinal chemical shift spectrum.

One step (the first step of this illustrative embodiment) is representedin the first column of FIG. 39 and is to form the reference spectra foreach phase group. Within this process, a first sub-step comprisesfrequency correcting each qualified frame by aligning the spectral peakto zero in the complex frequency domain, and to form the TD and FD meansfor each group. The next sub-step is to perform interpolated spectralpeak frequency shift estimation to determine the residual frequencyerror. This is done using cubic spline interpolation to increase theprecision of the estimate by a factor of 10, for example. This frequencyerror estimate is then used to correct the frequency of the TD mean foreach group. Because the preceding operations may have resulted in asmall phase shift, a final fine-grain phase adjustment is performed inthe TD. The final FD reference data is formed and both the TD and DFreferences for each phase group are available for spectralcross-correlation based frequency correction. Invoking the function“makePhGrpRef.m” in MATLAB, the commercially available processingutility by MathWorks, teaches how this is mechanized.

Another step, and which is step 2 of the particular illustrativeembodiment shown, is represented in the second column of FIG. 39,performs frequency correction using the spectral cross-correlationmethod on all frames. The first sub-step is to compute the complexspectral cross-correlation coefficient function between each frame andits corresponding group reference. This is done for a specified numberof lags, which may be typically for example about 200 for data fromabout a 3 Tesla system, sufficient to accommodate the largest frequency(for example about 50 Hz) shifts anticipated. Next, the correlation datais evaluated to locate the peak of the real part of the complexcorrelation function. Cubic spline interpolation is employed to increasethe precision by a factor of about 10 and the interpolated lag of thepeak is expressed as a frequency shift in Hz. Next, the frequency shiftestimation just determined is used in performing frequency correction onthe corresponding TD frame. The final step is to form the FD data foreach frame in each phase group. The MATLAB functions freqAdjSC.m,spectCorr.m, spectCorrinterp.m and cmnFreqAdjTD.m teach how theseoperations are mechanized.

Another step, which is step 3 in the illustrative embodiment shown, isrepresented in the third column of FIG. 39, and comprises forming thefinal spectrum by first combining the qualified frames by group in amanner to avoid introducing phase group cancellation bias due to theframe selection process. The mathematical means (or averages) ofqualified TD frames for each phase group are formed and then the mean ofthese means is computed. This provides the final (unless further stepsare desired) composite TD data from which the final chemical shiftspectrum is formed by applying appropriate apodization and transformingto the FD.

FIG. 39 shows a software flow diagram 500 of an automated method toperform frequency shift estimation and correction, and frame selectionin the context of phase group processing. It uses the spectralcross-correlation method of frequency shift estimation but could bebased on the spectral peak method of other approaches, including otherembodiments elsewhere herein described, just as well. The first function510 is to form reference spectra for each phase group. All the framesfrom the corresponding phase group will be aligned in frequency by thespectral cross-correlation method 520. The last step 530 associatesframe selection and phase group membership to form the mean of thefrequency-corrected and frame-selected frames for each phase group, andthen the mean of these means. This is performed in the time domain(“TD”) and then transformed into frequency domain (“FD”), withappropriate apodization applied to the FD used to form the finalchemical shift spectrum (or interim aggregate spectrum as the case maybe, to the extent other additional post-processing may be done, e.g.baseline correction, other filtering, curve fitting, etc.).

To form the reference spectra for each phase group 510, the first step511 according to the current illustrative embodiment is to frequencycorrect each “qualified frame” (frames which have passed frame editingas “selected” vs. “excluded” frames) using the peak grabbing & locationestimation method (such as for example, but not limited, according toother embodiments herein disclosed). That is used to form the TD and FDmeans (e.g. averages) for each phase group. The next step 512 is toperform interpolated spectral peak frequency shift estimation todetermine the residual frequency error. This is done using cubic splineinterpolation to increase the precision of the estimate by a factor of10. This frequency error estimate is then used to correct the frequencyof the TD mean for each group 513. Because the preceding operations mayhave resulted in a small phase shift, a final fine-grain phaseadjustment 514 is performed in the TD. The final FD reference data isformed 515 and both the TD and DF references for each phase group areavailable for spectral cross-correlation based frequency correction. TheMATLAB function makePhGrpRef.m teaches how this is mechanized.

Frequency correction using the spectral cross-correlation method 520 isnow performed on all frames. The first step is to compute the complexspectral cross-correlation coefficient function between each frame andits corresponding group reference 521. This is done for a specifiednumber of lags, typically 200 for data from a 3 Tesla system, sufficientto accommodate the largest frequency shifts anticipated. Next, thecorrelation data is evaluated to locate the peak of the real part of thecomplex correlation function 522. Cubic spline interpolation is employedto increase the precision by a factor of 10 and the interpolated lag ofthe peak is expressed as a frequency shift in Hz. Next, the frequencyshift estimation just determined is used in performing frequencycorrection on the corresponding TD frame 523. The final step is to formthe FD data for each frame in each phase group 524. The MATLAB functionsfreqAdjSC.m, spectCorr.m, spectCorrinterp.m and cmnFreqAdjTD.m teach howthese operations are mechanized.

The final spectrum is formed 530 by combing the qualified frames bygroup in a manner to avoid introducing phase group cancellation bias dueto the frame selection process. The means (averages) of qualified TDframes for each phase group are formed 531 and then the mean of thesemeans 532 is computed. This is the final composite TD data from whichthe final chemical shift spectrum is formed 533 by applying appropriateapodization and transforming to the FD.

FIGS. 40A-41D show raw spectral data to illustrate the benefits ofadvanced post processing to obtain usable results. Only phase correctionand apodization have been performed. In the following discussion, ‘fullspectrum” refers to an extended chemical shift range to include theresidual water signal. FIG. 40A shows the extended chemical shiftspectrum while “metabolite spectrum” shows an expanded plot of thechemical shift range in which important metabolite signals occur.

FIG. 40A shows the full spectrum of a single frame with an appropriatelevel of the residual water signal at chemical shift (CS) 4.7. Theprincipal metabolite of interest at CS 2.01 is hardly seen. Thisrepresents a generally good quality frame. FIG. 40B shows the metaboliterange for this same frame.

FIG. 40C shows the average of all 192 frames collected in this seriesand with only this very slight processing. Although the signal to noiseratio is greatly improved as expected by such averaging, the quality ofthe signal is greatly degraded by the inclusion of aberrant frames. Inthe metabolite range plot of the same spectrum shown in FIG. 40D, themetabolite peak which should be at CS 2.01 is broadened and shifted tothe left. This degradation reflects the need for frame quality screeningto exclude the aberrant frames and frequency correction to align thefrequency of the selected frames so they will combine constructively.

X. Baseline Correction

Baseline correction removes variations in the final spectral due tolarge macromolecules in the acquisition. Various approaches to baselinecorrection are contemplated as applicable with other embodiments of thisdisclosure, though they may not be expressly herein shown or described,as would be apparent to one of ordinary skill. One particular beneficialembodiment however is described by reference to the following baselinecorrection steps.

According to one step, the baseline is estimated using a “ranked orderfilter” (“ROF”) method. According to one further embodiment, the rankedorder filter is designed to find the local minimums of the spectra.However, other settings may be used to modify the portions or featuresof the spectrum tracked by or impacting the ROF results.

According to another step, outliers are removed from the ROF output. Theoutput of the ROF varies slowly and is influenced by large amplitudemetabolites. Adjacent points in the output of the ROF with relativelyhigh derivatives indicate large (and abrupt) changes. These will oftenrepresent real chemical peak regions, vs. baseline bias artifact.Accordingly, to distinguish baseline offset trends from real targetchemical peak offsets from baseline, and thus avoid “correcting” toflatten (and thus artificially reduce or lose) real peak regions to beredrawn along a baseline objective, points bounded by large amplitudemetabolites and points with high derivatives are considered outliers andare removed from the ROF output. This allows the next polynomial fitstep to track only on “non peak” offsets from a rolling baseline, andthus correct only for the baseline variances while preserving the peaksrising above it. A polynomial is fit to the ROF output without theoutliers. A polynomial is fit to the remaining points in the ROF. Inanother step, the polynomial is used to compute the baseline estimate.The baseline estimate is subtracted from the pre-baseline correctedspectrum to derive the final baseline corrected spectrum.

For further illustration, an example is provided by reference to FIGS.41A-4D as follows. FIG. 41A shows a post-processed spectrum from an MRSacquisition from a lumbar disc, baseline corrected using polynomial fitto an ROF incorporating all spectral data in the range (includingoutliers).

FIG. 41B shows certain baseline correction details overlaid to thespectrum, full spectrum range. More specifically, the ROF is shown belowthe spectrum in the metabolite range. The bold dotted line showspolynomial ft to the ROF without outliers, versus the narrower dashedline that shows the polynomial fit with outliers.

FIG. 41C shows these baseline correction details in a zoomed view morenarrowly drawn to a metabolite range of interest from the spectrum.

FIG. 41D shows a final spectra post baseline correction using thepolyfit without outliers contributing to the ROF used for the fit.

Additional examples that illustrate various of the foregoing aspects invarious combinations are provided as follows:

FIG. 42A shows a full range spectrum from an intravertebral disc MRSacquisition, post-processed according to a combination of multipleembodiments of the current disclosure.

FIG. 42B shows a time-intensity plot of the individual frames along theacquisition series providing the post-processed spectrum shown in FIG.42A, and shows prior to processing (top), with frame editing (middle),and after frequency correction of the retained qualified frames (bottom)(unqualified screened out frames replaced with dark horizontal lines inthe plot)

FIG. 42C shows a frame editing panel comprising a series of water signalanalyses comprising per-frame plots of “frequency estimate confidence”(confidence of peak location), frequency error, in-phase amplitude,FWHM, and finally qualified retained vs. excluded unqualified (red)frames at bottom.

FIG. 43A shows a similar spectral plot for the same acquisition shown inFIG. 42A, but bypassing frequency correction by setting the FWHMthreshold criteria below the FWHM of all frames (e.g. excluding them allthus bypasses frame editing and frequency correction all together).

FIG. 43B shows a time intensity plot for the same series of phases ofthe same acquisition shown in FIG. 43B.

FIG. 43C shows a frame editing panel comprising the water signalanalyses similar to that shown in FIG. 42C, and on the same acquisitionfeatured in FIG. 42C, but related to the different frame editing andfrequency correction thresholds and related processing featured in FIGS.43A-B.

While certain embodiments of the disclosure have been described, theseembodiments have been presented by way of example only, and are notintended to limit the scope of the broader aspects of the disclosure.Indeed, the novel methods, systems, and devices described herein may beembodied in a variety of other forms. For example, embodiments of oneillustrated or described MRS system post-processing component may becombined with embodiments of one or more illustrated or described MRSsystem processing components. Moreover, the MRS system componentsdescribed herein, e.g. pulse sequence, signal processor, or diagnosticprocessor, may be deployed for particular beneficial use forintervertebral discs, or utilized for other purposes. For example, anMRS system (or component sequence, signal processor, or diagnosticprocessor useful therewith or therein), may be configured and used inmanners consistent with one or more broad aspects of this disclosure fordiagnosing other tissue environments or conditions than pain within anintervertebral disc. Or, such may be usefully employed for diagnosingpain or other tissue environments or conditions in other regions ofinterest within the body. Such further applications are consideredwithin the broad scope of disclosure contemplated herein, with orwithout further modifications, omissions, or additions that may be madeby one of ordinary skill for a particular purpose. Furthermore, variousomissions, substitutions and changes in the form of the methods,systems, and devices described herein may be made without departing fromthe spirit of the disclosure. Components and elements may be altered,added, removed, or rearranged. Additionally, processing steps may bealtered, added, removed, or reordered. While certain embodiments havebeen explicitly described, other embodiments will also be apparent tothose of ordinary skill in the art based on this disclosure.

In one further embodiment contemplated, a computing system, comprisingone or more microprocessors receiving at least one signal responsive todata collected in an MR scanner, is configured to implement a magneticresonance spectroscopy (MRS) processing system configured to process arepetitive multi-frame MRS spectral acquisition series generated andacquired for a voxel principally located within an intervertebral discvia an MRS pulse sequence, and acquired at multiple parallel acquisitionchannels of a multi-coil spine detector assembly, in order to providediagnostic information associated with the disc. According to one stillfurther embodiment, this comprises: an MRS signal processor comprising achannel selector, a phase shift corrector, a frequency shift corrector,a frame editor, and a channel combiner. The system of this embodiment isconfigured to receive and process the MRS spectral acquisition seriesfor the disc and to generate a processed MRS spectrum for the serieswith sufficient signal-to-noise ratio (SNR) to acquire informationassociated with identifiable features along MRS spectral regionsassociated with unique chemical constituents in the disc. An MRSdiagnostic processor is configured to extract data from identifiablechemical regions in the processed MRS spectrum in a manner that providesdiagnostic information for diagnosing a medical condition or chemicalenvironment associated with the disc.

In another embodiment, a non-transitory physical computer readablemedium encoded with a computer executable computer program or code isoperable by a processor to run the computer program or code to cause acomputing system to implement a magnetic resonance spectroscopy (MRS)processing system configured to process a repetitive multi-frame MRSspectral acquisition series generated and acquired for a voxelprescribed to correspond with a tissue region of interest in the body ofa patient. The computer program or code comprises one or more (invarious combinations) of the processing embodiments disclosed hereunder.According to a further embodiment, the voxel is principally locatedwithin an intervertebral disc via an MRS pulse sequence, and acquired atmultiple parallel acquisition channels of a multi-coil spine detectorassembly, in order to provide diagnostic information associated with thedisc. According to still another further embodiment, the computerprogram or code comprises an MRS signal processor that further comprisesone or more (in combination) of the following: a channel selector, aphase shift corrector, a frequency shift corrector, a frame editor, anda channel combiner. This is configured to receive and process the MRSspectral acquisition series for the disc and to generate a processed MRSspectrum for the series with sufficient signal-to-noise ratio (SNR) toacquire information associated with identifiable features along MRSspectral regions associated with unique chemical constituents in thedisc. An MRS diagnostic processor is provided under still anotherfurther embodiment, and is configured to extract data from identifiablechemical regions in the processed MRS spectrum in a manner that providesdiagnostic information for diagnosing a medical condition or chemicalenvironment associated with the disc.

In another further embodiment, a magnetic resonance spectroscopy (MRS)processing method is used for processing a repetitive multi-frame MRSspectral acquisition series generated and acquired for a voxelprincipally located within an intervertebral disc via an MRS pulsesequence, and acquired at multiple acquisition channels of a multi-coilspine detector assembly, and for providing diagnostic informationassociated with the disc. According to certain further embodiments, theMRS processing method comprises one or more (in combination) of theother processing component embodiments disclosed hereunder. The methodaccording to another particular embodiment comprises: receiving the MRSspectral acquisition series from the multiple acquisition channels; andsignal processing the MRS acquisition series. The signal processingcomponent of this method comprises: selecting one or more channels amongthe multiple channels based upon a predetermined criteria, estimatingand correcting phase shift error among multiple frames within the seriesof a channel acquisition, estimating and correcting a frequency shifterror between multiple frames within the series of the channelacquisition, flagging and editing out sub-optimal frames from the seriesbased upon a predetermined criteria frame selection criteria, combiningselected and corrected channels for a combined average processed MRSspectrum. Another method embodiment comprises adjusting spectral signalupfield of water signal by subtracting artifact derived from downfieldsignal components generally mirrored opposite about a water peak centerline. Still further method embodiments comprise diagnosticallyprocessing the processed MRS spectrum by extracting data fromidentifiable chemical regions in the processed MRS spectrum andprocessing the extracted data in a manner that provides MRS-baseddiagnostic information for diagnosing a medical condition or chemicalenvironment associated with the disc.

It is further contemplated that other processing approaches known in theart may be also combined with the aspects, modes, and embodimentsexpressly shown and described hereunder, such as for example eddycurrent correction, other apodization techniques etc. Moreover, othertypes of MRS acquisitions, such as multi-voxel for example, may be used.Still further, the processors and related methods disclosed may beembedded into the controller or processing environments of the MRsystems themselves used in the acquisitions, or may be stored, operated,and conducted remotely. For example, such may be stored and/or operatedvia a separate server, either co-located or remote from the MRacquisition system, such as via data received electronically via theinternet, PACS, FTP uploading, or tangible physical media such as mobilehard drives or CDs/DVDs. It is further contemplated that notebookcomputers, PDAs, mobile phones, etc. may also provide at least in partthe processing, communicating, and/or display functions according to thepresent aspects, modes, and embodiments of this disclosure.

TABLE 1 Range of Optimization Thresholds. Number of Parameter Range ofValues Values Confidence Threshold 0.5:0.1:0.7 3 FWHM Threshold0.15:0.05:0.35 5 Frequency Error 0.15:0.05:0.30 4

What is claimed is:
 1. A method for processing a set of multipleserially acquired magnetic resonance spectroscopy (MRS) free inductiondecay (FID) frames from a multi-frame MRS acquisition series from aregion of interest (ROI) in a subject, and for providing apost-processed MRS spectrum, comprising a combination of at least two ofthe following: performing the step on the set at each of the multipleparameter values and providing multiple processed results, respectively,comprising multiple respective feature values, respectively, comparingthe feature values of the processed results, determining a chosenfeature value among the multiple feature values based upon thecomparison and corresponding to a chosen parameter value, setting theparameter to the chosen parameter value for performing the processingstep to provide the post-processed MRS spectrum; providing a firstpost-processed MRS spectrum, identifying a feature of a downfield regionof the spectrum that is located downfield along a larger part permillion (ppm) range than a water peak location, comparing the featurewith a reference, performing a processing operation based on thecomparison and related to the post-processed MRS spectrum; determining afrequency shift error of an FID frame by performing spectralcross-correlation between the FID frame and a reference frame; adjustinga phase of an FID frame by a phase shift configured to reduce an extentof the absorption spectrum having a negative value; estimating a linewidth between two opposite walls of a water peak region of an FID frameof the series, comparing the estimated line width against a thresholdvalue, and qualifying and including, or disqualifying and excluding, theFID frame from the set for further processing to provide thepost-processed MRS spectrum based at least in part upon the comparison;associating first and second groups of FID frames within the set withfirst and second different relative phases, respectively, of a phasestep cycle along the MRS acquisition series, processing the first andsecond groups separately, averaging the first and second groupsseparately to form first and second average interim spectra, andaveraging the first and second average interim spectra; and generating afirst baseline estimate using a rank order filter (ROF) on a range ofinterest of an interim MRS spectrum partially post-processed from theset, fitting a polynomial to at least a portion of the first baselineestimate to generate a baseline curve, and subtracting the baselinecurve from the interim MRS spectrum to provide the post-processed MRSspectrum.
 2. One or more non-transitory computer readable mediacomprising computer instructions configured to cause one or morecomputer processors to perform actions comprising the method of claim 1.